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SECTION  I 


INTRODUCTION 


The  past  decade  has  seen  an  order-of-magnitude  increase  in  reli¬ 
ance  on  fixed-  and  moving-base  simulation  for  research  and  training  in 
all  types  of  aircraft.  The  reasons  for  this  increase  are  numerous, 
ranging  from  the  development  of  vehicles  which  must  he  completely 
design-validated  before  flight  test  (e.g.,  the  Space  Shuttle)  to  the 
increased  use  of  simulation  for  training  purposes  (for  reasons  of 
safety,  repeatability,  and  reduction  of  fuel  costs).  As  a  result,  simu¬ 
lation  facilities  are  in  widespread  use  throughout  the  Air  Force  as  well 
as  other  government  agencies  and  private  industry. 

Accompanying  this  increased  use  of  simulation  has  been  a  dramatic 
increase  in  simulator  complexity.  This  particular  development  manifests 
itself  most  noticeably  as  digital  computers  replace  analog  equipment. 
While  digital  implementations  provide  better  reliability  and  mainten¬ 
ance,  increased  static  accuracy,  and  greater  flexibility,  they  also 
introduce  an  array  of  simulation  artifacts  heretofore  unconsidered. 
These  anomalies  impact  the  new  generation  of  simulators  in  two  areas: 

•  Hardware/software  procurement.  Techniques  used 
to  specify  analog  systems  will  no  longer  suffice 
in  a  digital  environment.  Critical  concerns  such 
as  frame  time,  word  length,  integration  algo¬ 
rithms,  data  skewness,  order  of  subroutine  call, 
etc.,  are  presently  the  heuristic  choice  of  the 
contractor,  since  the  contracting  agency  has  no 
analytical  means  by  which  to  specify  these  items. 

•  Research/training  on  existing  equipment.  Often 
digital  simulation  artifacts  creep  into  a  simula¬ 
tion  facility  as  a  result  of  upgrading  existing 
equipment  with  digital  computers.  As  a  result, 
experiments  and  training  sessions  may  be  contam¬ 
inated  with  extraneous  simulation  errors.  These 
errors  are  difficult  to  detect  and  assess  without 
the  aid  of  analytical  tools. 
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The  anomalies  them  selves  are  many  but  can  be  roughly  described  as 
frequency  aliasing  effects  (another  term  is  folded  power).  They  arise 
for  a  variety  of  reasons: 

1)  Two  or  more  computers  required  in  a  large  simula¬ 
tion,  each  working  in  its  own  frame  time  (the 
so-called  Independent  processor  problem 

2)  Serial  processing  (calling)  of  subroutines.  The 
first  subroutine  called  may  work  with  different 
input  data  than  those  called  later  (skewed  data, 

"stale"  input  data). 

3)  Throughput  delay  factors. 

4)  Staircasing  (zero-order-hold  effects)  when  the 
digital  computer  output  is  used  to  drive  the 
actuators  of  motion-base  cabins. 

A  set  of  recently  developed  concepts  provides  the  basis  for  iden¬ 
tifying  potentially  critical  simulation  anomalies  at  the  design  stage  in 
an  organized  and  rational  way.  Moreover,  a  method  exists  which  can  be 
used  to  implement  given  computer  code  (integration  algorithms)  into  a 
multi-rate,  time  delayed,  skewed  data  analytical  model,  and  predict  the 
impact  of  these  digital  effects  on  the  proposed  simulation.  To  date 
these  new  methods  have  been  demonstrated  on  low-order  models  with  only 
limited  demonstration  of  the  theory.  It  is  necessary  to  demonstrate  all 
features  of  the  new  theory  in  a  joint  fashion  on  higher-order  problems- 

The  methodology  embodies  three  concepts.  The  first,  which  will  be 
illustrated  in  Section  II,  is  the  continuous  frequency  response  of  a 
digitally  controlled  system.  Using  techniques  developed  in  References  1 
and  2,  it  is  possible  to  compute  the  group  of  N  sinusoids  which  fit  the 
response  of  a  digitally  controlled  system  not  only  at  the  sample  points 
but  at  the  (N  -  1)  inter-sample  points  as  well.  In  the  limit,  as  N 
approaches  infinity,  one  obtains  the  "continuous"  frequency  response  of 
a  digitally  controlled  system.  This  section  will  clarify  the  term  "fre¬ 
quency  response"  in  the  content  of  this  report;  it  will  be  demonstrated 
that  there  is  a  truly  significant  difference  between  the  frequency 
response  of  a  digitally  controlled  continuous  system  and  the  discrete 
spectrum  of  sampled  data  control  theory.  It  is  convenient  to  review  the 


theory  of  the  single-rate  frequency  response  since  it  affords  the  oppor¬ 
tunity  to  present  several  clarifying  comments  on  the  frequency  response 
concept  suggested  by  various  readers  of  Reference  1  .  Thus,  Section  II 
is  also  an  "update"  of  the  frequency  response  concept  for  digitally  con¬ 
trolled  systems  and  provides  a  link  with  newer  developments. 

The  second  component  of  this  simulation  analysis  method  is  vector 
switch  decomposition,  a  technique  for  analyzing  simulations  with  two  or 
more  sample  rates,  data  skewness,  throughput  delays,  etc.  As  described 
in  Section  III,  this  technique  is  conceptually  quite  simple.  All  samp¬ 
lers  in  a  given  system  are  replaced  by  equivalent  samplers  whose  periods 
are  the  least  common  sampling  period,  and  the  appropriate  time  delay 
vectors.  One  practical  problem  with  this  method  is  the  accompanying 
increase  in  dimensionality  for  the  "decomposed"  vector.  Vector  switch 
decomposition  is  not  a  particularly  workable  tool  for  pencil  and  paper 
design,  although  it  is  quite  amenable  to  computerization,  since  the 
matrix  manipulations  are  routine. 

We  then  develop  from  switch  decomposition  an  algebra  that  circum¬ 
vents,  for  a  limited  class  of  problems,  the  dimensional  complexities 
introduced  by  the  decomposition  itself.  In  effect,  a  scalar  problem 
will  remain,  in  the  framework  of  this  algebra,  a  completely  scalar  prob¬ 
lem-  This  "limited  class  of  problems"  is  important  since  it  encompasses 
the  open-loop  analysis  of  particular  elements  of  a  simulation  as  well  as 
closed-loop  multi-rate  systems  wherein  the  sample  rates  are  related  as 
powers  of  2.  This  latter  case  covers  both  the  Space  Shuttle  (25,  50, 
100  Hz)  and  the  F-18  (20,  40,  and  80  Hz)  digital  control  systems. 

An  important  computational  aspect  of  this  scalar  algebra  is  dis¬ 
cussed  extensively  in  Section  IV.  Specifically,  an  algorithm  is  des¬ 
cribed  for  transforming  from  a  T/N  time  frame  to  a  T/M  time  frame,  N  and 
M  being  arbitrary  but  rational.  For  example,  sets  of  (M,  N)  such  as 


*In  particular,  the  comments  of  Dr.  Hsi-Han  Yeh  of  the  University  of 
Kentucky  provided  an  excellent  interpretation  of  the  limiting  results  as 
N  ♦  infinity. 
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(3,  2),  (2,  3),  (7.8,  1),  etc.,  are  permissible.  The  algorithm  is  of 

particular  importance  as  it  permits  the  user  to  circumvent  certain 
dimensionality  constraints  introduced  by  multi-rate  methods  based  on 
residue  theory.  The  nature  of  these  constraints  will  be  made  clear  in 
the  case  studies  of  Sections  VII  through  IX. 

In  Section  V,  the  application  of  the  scalar  algebraic  approach  to 
multi-rate  frequency  response  is  compared  with  the  more  general  switch 
decomposition  approach  of  Reference  2. 

The  third  component  of  the  simulation  analysis  method  was  developed 
in  Reference  2.  It  is  primarily  a  technique  to  incorporate,  within  the 
switch  decomposition  framework,  a  specified  computer  code.  For  example, 
there  are  a  variety  of  methods  available  for  modeling  a  low-pass  filter 
section  on  a  digital  computer.  One  can  use  the  Tustin  transform,  a  rec¬ 
tangular  integration  algorithm,  Adams-Bashforth,  etc.  Clearly,  in 
analyzing  a  given  simulation  there  must  be  a  capability  for  incorporat¬ 
ing  the  given  difference  equations  into  the  analysis,  without  any  par¬ 
ticular  regard  (or  prejudice)  as  to  what  the  originator  intended  the 
code  to  represent.  It  is  the  task  of  the  analysis  to  show  the  origin¬ 
ator  how  successful  his  digital  model  is  in  representing  the  principal 
features  of  the  continuous  system.  Facets  of  this  problem  are  discussed 
in  Section  VI. 

In  Section  VII,  a  first  effort  is  made  to  bring  all  of  the  key  ele¬ 
ments  together  into  a  two-rate  simulation  case  study.  The  ratio  between 
the  two  frame  times  forces  the  use  of  the  switch  decomposition  format 
and  furnishes  insight  into  the  dimensionality  problems  encountered. 
This  case  study  also  demonstrates  how  a  two  rate  format  can  introduce 
unintended  lightly  damped  modes  into  a  simulation. 

In  Section  VIII,  another  case  study  is  described  which  does  not 
require  the  use  of  switch  decomposition.  The  primary  purpose  of  this 
study  is  to  gain  insight  into  the  multi-rate  frequency  response.  Sup¬ 
pose  the  output  is  sampled  in  a  T/3  time  frame  but  other  rates  are 
Involved  in  preceding  portions  of  the  system.  How  many  sine  waves  are 
required  to  exactly  match  the  steady  state  T/3  output  samples? 


4 


In  Section  IX  the  analysis  of  an  existing  three  rate  simulation  is 
attempted.  It  was  in  this  high-order  simulation  study,  where  the  ratios 
between  frame  times  were  very  large,  that  shortcomings  in  our  computa¬ 
tional  tools  proved  to  be  more  critical  than  previously  suspected.  For 
instance,  in  moving  from  one  time  frame  to  another,  a  fourteenth-order 
system  becomes  a  system  of  order  1121  The  "invention"  of  the  algorithm 
of  Section  IV  was  a  direct  result  of  these  difficulties.  Subsequently, 
we  were  able  to  proceed  through  the  case  study  with  relative  ease  and 
achieve  important  results.  Thus  Section  IX,  as  it  now  stands,  will  give 
little  insight  into  the  numerical  and  dimensional  difficulties  first 
encountered . 

The  multi-rate  simulation  studies  of  Sections  VII  through  IX  have  a 
primary  goal  of  providing  insight  into  the  spectral  characteristics  of 
the  output  steady-state  wave form  in  terms  of  the  number  of  sinusoids 
required  to  match  the  output  (sampled)  data  points.  A  secondary  objec¬ 
tive  is  to  call  attention  to  some  observed  anomalies  (such  as  extraneous 
lightly  damped  modes  introduced  by  the  multi-rate  structure)  which  can 
be  identified  and  quantified  using  the  anlaytical  tools  described  in  the 
earlier  sections  of  the  report. 

Section  X  treats  a  single-rate  case  study  which  investigates  the 
effective  and  unintended  filtering  introduced  when  subroutines  are 
called  in  a  serial  manner.  Specifically,  it  is  shown  that  the  z-domain 
analytical  model  of  the  computer  code  used  in  the  implementation  of  a 
washout  filter  for  a  large  moving-base  simulator  is  twice  the  order  of 
the  intended  s-domain  transfer  function. 

The  report  is  a  lengthy  one,  attributable  in  a  large  part  to  a 
desire  to  pull  together  under  one  cover  the  key  elements  of  References  1 
and  2  which  have  application  to  simulation  error  analysis.  Those 
readers  familiar  with  the  concepts  of  References  1  and  2  are  therefore 
in  a  position  to  selectively  read  the  present  report. 
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SECTION  II 


DISTINCTIONS  IN  "FREQUENCY  RESPONSE" 


A.  INTRODUCTION 

The  term  "frequency  response"  for  discrete  systems  refers  to  the 
process  of  selecting  the  magnitude  and  phase  of  a  single  frequency  sine 
wave  to  fit  sample  points  at  the  sampling  instants*  In  contrast,  the 
concept  of  "frequency  response"  for  the  continuous  variables  of  a  dis¬ 
cretely  (digitally)  controlled  system  defines  the  infinite  set  of  sine 
waves  (In  terms  of  a  fundamental  and  aliased  frequencies)  which  add 
together  to  exactly  reproduce  the  continuous  output  steady  state  wave¬ 
form. 

A  simple  but  illuminating  example,  defined  by  Figure  1,  illustrates 
the  distinction.  Suppose  the  system  is  forced  by  a  step  input  and  the 
continuous  output  C  is  recorded  —  both  as  a  continuous  analog  quantity 
and  also  sampled  at  a  rate  of  i/T  Hz.  Furthermore,  let  the  open-loop 
plant,  10/(s  +  10),  be  subject  to  two  different  control  laws,  one  imple¬ 
mented  with  a  zero-order  hold  (ZOH),  the  other  with  the  "slewer"  data 
hold.  The  control  laws  were  synthesized  with  the  objective  of  forcing 
the  output  to  be  the  same  at  the  sampling  instants,  regardless  of  the 
control  law/coupler  being  used.  These  responses  are  shown  in  Figure  2. 
Note  the  smoothness  of  the  slewer-controlled  response  and  the  roughness 
of  the  ZOH  response.  However,  an  observer  who  is  shown  only  the  1/T  Hz 
sampled  output  would  be  unable  to  detect  any  differences  in  the  tran¬ 
sient  responses,  even  though  the  continuous  responses  are  very  differ¬ 
ent. 

Next  force  each  system  with  a  sine  wave  and  record  both  the  continu¬ 
ous  and  discrete  output  waveforms.  The  observer  who  is  only  interested 
in  matching  the  sample  points  uses  the  discrete  frequency  response  and 
picks  the  magnitude  and  phase  from  a  Bode  plot,  such  as  the  one  labeled 
"discrete,"  in  Figure  2b.  Again,  this  observer  is  unaware  of  any  dis¬ 
tinction  between  the  two  systems  (ZOH  or  slewer)  as  the  same  sine  wave 
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fits  the  sampled  points  of  either  output  waveform.  The  "continuous" 
frequency  response  magnitude  plots  for  the  ZOH  and  slewer  designs  are 
also  shown  In  Figure  2b.  Here  the  observer  must  know  how  to  reproduce 
the  output  waveform  as  the  sum  of  a  fundamental  and  Its  aliased  fre¬ 
quencies  —  a  point  discussed  next.  We  merely  observe  that  the  "con¬ 
tinuous"  Bode  plots  show  a  truly  significant  difference  between  the  two 
systems  —  differences  which  the  discrete  frequency  response  of  classi¬ 
cal  sample  data  control  theory  is  Incapable  of  detecting. 

In  the  analysis  of  multi-rate  simulations,  one  is  often  more  inter¬ 
ested  in  the  finite  set  of  sine  waves  that  fits  an  output  sequence 
sampled  In,  for  example,  a  T/M  time  frame  when  the  input  Is  sampled  In  a 
T/N  time  frame.  It  Is  therefore  appropriate  to  review  the  "finite  N 
case"  and  the  subsequent  extension  to  the  continuous  frequency  response 
of  a  discretely  excited  system. 

B.  FREQUENCY  RESPONSE  OF  A  SAMPLED  SYSTEM 

When  a  continuous,  stable,  linear  system  is  excited  by  a  sine  wave, 
the  steady-state  waveform  Is  comprised  of  a  single  wave  at  the  same  fre¬ 
quency  as  the  Input.  It  differs  from  the  Input  wave  only  by  a  phase 
angle  and  a  magnitude  factor.  Moreover,  it  Is  unnecessary  to  compute 
the  actual  transient  response  of  the  system  when  the  behavior  for  large 
values  of  time  is  of  Interest,  since  both  the  magnitude  factor  and  phase 
angle  can  be  read  from  a  Bode  plot. 

The  analysis  of  a  continuous  system  which  is  discretely  controlled 
is  more  complex.  Given  that  the  system  is  stable,  the  continuous  output 
waveform  will  contain  more  than  just  a  wave  at  the  fundamental  fre¬ 
quency.  It  will  consist  of  the  fundamental  and  all  of  its  aliases. 
Thus,  if  the  system  is  forced  with  l  sin  bt,  the  output  will  contain 
terms  at  frequencies  b,  (b  +  (2x/T)J,  [b  -  (2*/T)J,  ....  The  relative 
amplitudes  and  phase  angles  will  depend  on  the  data  hold  employed  as 
well  as  the  system  transfer  function.  Nevertheless,  given  a  data  hold 
and  transfer  function,  the  magnitude  and  phase  angle  of  each  component 
can  be  read  from  a  particular  Bode  plot.  This  concept  of  frequency 
response  Is  more  comprehensive  than  the  traditional  concept  of  the 
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"sampled  spectrum,"  which  is  limited  to  determining  the  single  sinusoid 
that  fits  the  system  output  samples  at  the  sampling  instants. 

In  the  subsections  to  follow  the  Bode  plot  concept  for  a  continuous 
system  is  reviewed  and  then  extended  to  the  frequency  response  of  a  dis¬ 
cretely  excited  open-loop  system. 

C.  CONTINUOUS  SYSTEM  BODE  PLOTS 

It  will  be  helpful  to  first  review  the  Bode  plot  concept  for  contin¬ 
uous  sytems.  Let  R  in  Figure  3  be  a  unit  input  sine  wave  with  frequency 
u)Q  rad/sec.  The  output,  in  the  frequency  domain,  is: 

U) 

C(s)  =  G(s)R(s)  =  G  (s)  (1) 

sz  +  u>g 

Equation  1  can  be  expanded  in  partial  fractions  as: 


C  (s) 


Aai0  Bs 

s  ^  s^  +  u)£ 


+ 


Terms  associated 
with  characteristic 
polynomial  of  G(s) 


(2) 


Given  that  all  poles  of  G(s)  are  in  the  left  half  plane,  the  bracketed 
term  in  Equation  1  represents  time  functions  that  vanish  as  t  -*■ 
Thus,  the  steady-state  behavior  is  completely  determined  by  the  partial 
fraction  coefficients  A  and  B,  and  once  they  are  known  the  steady-state 
time  response  can  be  written  directly  as: 


Figure  3.  Continuous  System 
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C(t) 


11m 


A  sin  0)ot  +  B  cos  w0t 


I  £  *00 


A 


A  +  B  sin  (u)_t  +  $) 


(3) 


where  $  *  tan  *  (B/A) • 

The  details  of  solving  for  A  and  B  show  the  relationship  between  the 

Bode  plot  and  the  steady-state  waveform.  To  solve  for  A  and  B,  multiply 
2  2 

Equation  2  by  [s  +  u>0]  and  evaluate  the  result  for  s  -  jw0S 

G(s)w0|  .  -  (Au)0+Bs)|  + 

°|s-jo)0  °  |s-jtu0 

(4) 

or 


Terms  associated 
with  characteristic 
polynomial  of  G(s)  I 


s 2+u2) | 


S“jwQ 


A  +  JB 


/A2  +  B2  tan_1  B/A 


(5) 


To  summarize,  a  sinusoidal  input  at  frequency  u>0  produces  a  steady-state 
waveform  of  the  same  frequency.  It  differs  from  the  input  only  by  a 
magnitude  factor  and  a  phase  shift.  Both  the  magnitude  factor  and  phase 
shift  for  any  given  input  frequency,  u>o,  can  be  read  directly  from  the 
Bode  plot  for  G(ju)).  That  is,  for  any  given  input  frequency 

A  +  IB  -  G(s)  I  ,  (6) 

This  "frequency  response"  viewpoint  is  expanded  to  encompass  discretely 
excited  continuous  systems. 

D.  MATHEMATICAL  PRELIMINARIES 

Let  R  be  a  sinusoid  of  unit  amplitude  with  frequency  b  rad/sec 
(R  -  1  sin  bt).  If  R  is  sampled  at  1/T  Hz  and  then  described  in  terms 
of  an  N/T  Hz  model,  we  obtain 
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**> 


rt 


_ sin  bT _ 

z2N  -  (2  cos  bT)zN  +  1 


2  -  eaT/N 


(7) 


T 

where  the  superscript  [*]  Is  used  to  denote  the  frametime  of  the  sampl¬ 
ing  operator*  For  later  use,  it  is  necessary  to  find  the  2N  factors  of 
the  denominator  of  Equation  7  in  a  form  that  will  permit  a  partial  frac¬ 
tion  expansion  containing  terms  for  which  corresponding  time  functions 
are  known.  For  example,  if  f(t)  «  sin  bt,  then 


fT(t)  *  [sin  bt]T  =>  F(z)  =  ~x - ~S~n— ^ - 

zz  -  (2  cos  bT)z  +  1 


(8) 


but  we  do  not  know  the  time  function  corresponding  to 


F(z) 


_ _ z  8  in  bT _ 

z2  +  (2  cos  bT)z  +  1 


(9) 


which  will  also  occur  among  the  N  factors  of  the  denominator  of  Equa¬ 
tion  7. 

N 

This  problem  can  be  examined  in  more  detail  by  letting  z  =  X  in 
Equation  7  and  solving  for  the  roots  using  exponential  notation: 

z2N  -  (2  cos  bT)zN  +1  -  X2  -  2(cos  bT)X  +  1  (10) 


In  turn. 


X2  -  (2  cos  bT)X  +  1 


(X  -  cos  bT)2  +  (sin  bT)2 


(H) 


The  N  roots  of  Equation  11  can  now  be  expressed  as 


(X  -  cos  bT)2  -  -(sin  bT)2 
X  -  cos  bT  -  +j  sin  bT 

X  ■  cos  bT  +  j  sin  bT  ” 

Replacing  X  by  z  gives  one  of  the  many  ways  of  describing  the  roots  of 
Equation  7: 

z  -  e+JO>T/N)  _  e+j{(bT/N)+(2irn/N)]  ,  [(bT/N)  +  (2*n/N)]  (12) 

In  a  purely  formal  sense,  the  n  in  Equation  12  can  take  on  both 
positive  and  negative  integer  values.  The  preferred  format  for  defining 
the  roots  of  Equation  12  is: 


As  we  have  said,  both  positive  and  negative  values  of  n  are  permitted. 
For  example,  if  N  -  3  there  are  three  complex  conjugate  roots  pertaining 
to  the  frequencies 


b,  b  +  2*/T  ,  b  +  4n/T 


However,  the  values 


b,  b  -  (2x/T)  ,  b  +  (2x/T) 

are  equally  permissible.  For  the  finite  N  case  many  readers  will  prefer 
the  description  in  terms  of  the  input  frequency  and  the  positive  ali¬ 
ases,  avoiding  a  description  that  contains  negative  frequencies. 
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S.  OPEN-LOOP  FREQUENCY  RESPONSE  —  FINITE  N 


Consider  the  system  of  Figure  4  where  G(s)  represents  an  arbitrary 
transfer  function  and  M  represents  an  arbitrary  data  hold.  Suppose  R  is 
a  unit  amplitude  sine  wave  and  the  output  is  sampled  in  a  T/N  frame 
time : 


(GM)T/Nrt  .  (GM)T/N 


sin  bT 


Z2N  _  ( 2  cos  bT)z^  +  1 


The  superscript  is  used  to  call  out  the  defining  time  frame  and  use  is 
made  of  the  identity  [a  B^-j  «  a^^rT  (Reference  1). 

Expand  the  right-hand  side  of  Equation  13  in  partial  fractions: 


N-l  Anz  sin  Mnd/N)  +  Bnz[z  -  cos  u>n(T/N)] 
n=0  z2  -  [2  cos  uin (T/N)  ]  z  +  1 


+  [Terms  due  to  modes  of  (GM)T/NJ  (14) 

T /N 

Assume  that  responses  in  the  modes  of  (GM)  approach  zero  as  t  ■+■  00 , 
i.e.,  that  all  modes  are  stable.  In  Equation  14, 


*9 


m 


For  the  present,  we  assume  that  b  <.  2»T.  The  steady-state  waveform,  at 
the  sampling  Instants,  can  be  written  as: 


[c(t)]T 


T 

N-l 

T  (A„  sin  +  Bn  cos  “>nt) 

.  n“0 


(15) 


To  solve  for  and  Bn,  multiply  each  side  of  Equation  14  by 
[z2  -  [2  cos  <*  (T/N)J z  +  l]  ,  0  <  k  <  (H  -  1) 

and  evaluate  for  z  -  14<VT/N>'  The  onl?  term  that  030  ^  ^ 

right-hand  side  occurs  when  n  -  k  (the  k  notation  can  then  be  changed  to 

n) .  To  illustrate. 


(GM)T/N  - zLain.tl.- - -  [z2  _  2  cos  u>k(T/N)z  +  D  *-U<«k(T/N) 

z2N  .  (2  coa  bT)zN  +1  '  * 


[Anz  sin  ^(T/N))  +  Bnz[z  -  cos 

*2  -  [2  cos  ujjjd/N)]  z  +  1 


(2  cos  u^CT/tOJz  +  1]  t.Wuk(T/B) 


(16) 


For  any  n  ^  k,  the  right-hand  side  of  Equation  16  is  identically  zero 
since 

z2  -  [2  cos  a>k(T/N)]z  +  1  -  [z  -  cos  ^(T/N)]2  +  lain  ^(T/N)]2  (17) 

vanishes  when 

z  -  14u)k(T/N)  -  cos  wk(T/N)  +  j  sin  wk(T/N)  (18) 
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Specifically,  we  obtain 


r\  r% 

[cos  uk(T/N)  +  j  sin  uik(T/N)  -  cos  uik(T/N)]  +  [sin  u>k(T/N)j  i  0 

(19) 

For  n  “  k,  the  cancellation  of  the  common  factor  guarantees  the  survival 
of  an  n  «  k  term.  Factoring  out  a  common  2  gives 

z[Ak  sin  oon(T/N)  +  Bk(cos  wn(T/N)  -  cos  u»n(T/N)  +  j  sin  uin(T/N)]j 

=  (Ak  +  jBk)z  sin  “k(T/N)|z.Uuk(T/N) 


Therefore,  Equation  16  becomes 


(20) 


z  sin  ->k(T/N)  (Ak  +  JBk) 


(Q/H)T/s(ts  sin  bt)(z2  -  [2  cos  u>k (T/ff Jlc  +  l) 
z2N  -  12  cos  bT)zN  +  1 


z-l»wk (T/N; 


(21) 


At  this  point,  let  k  revert  to  n. 


(CM)T/N j 


.  zN-1  sin  hT 
■Uwn(T/N)  sin  wn(T/N) 


z2  -  12  COS  UndVNHz  +  l) 

z2N  -  (2  cos  bI)zN  +  1  c-l^nCr/N) 


(22) 


The  last  term  on  the  right-hand  side  of  Equation  22  is  indeterminate 
(0/0)  when  z  -  Hu>n(T/N).  Therefore,  apply  L'Hopital's  rule  once  and 
obtain 


*1.  +  JBn 


(CM)T/N 


*-14Un(T/N) 


-  COS  inn(T/N)l  Bln  bT  I 
2Nzn~1(zn  -  cos  bT)  •*» 


(23) 


4^- 
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+  J“n 


-  i  (CM)  T/N 

N 


z-U-pd/N) 


!  cos  un  (T/N)  +  J  Sin  uin(T/N)  -  eoa  u>„(T/N)]  Bln  bT_ 
(cos  ^nT  +  ~j  sln  u.nT  -  cos  bTj  8ln  “‘n<T/N) 


(24) 


A  direct  substitution  for  .  b  +  (2*n/T)  quickly  shows  that  the  last 
part  of  the  product  in  Equation  24  Is  unity*  Therefore, 


+  JBn 


±  (GM)T/H 
N 


*  ST/N 
z-e 

z-14u.n(T/N) 


(25) 


The  superscript  notation  In  Equation  25  is  for  the  purpose  of  calling 
out  the  definition  of  z  being  used  in  the  evaluation. 

To  review  the  situation,  the  system  is  forced  with  the  sin  bt*  The 
steady-state  output  waveform,  sampled  with  a  T/N  frametime,  has  the  form 


CT/N(t) 


N-l 


E  (An  sin  uint  +  Bn  cos  u>nt) 


T/N 


where 


b  + 


2*n 


0,  1,  2, 


N-l 


(26) 


(27) 


As  an  alternative,  one  may  use 


CT/N  - 


(N-D/2 

E  (An  sin  u»nt  +  ^  cos  u>nt) 
In— (N-l)/2 


T/N 


,  N  odd  (28) 


or 


CT/N 


(N  /  2 )  —  1 

E  (An  sin  u)nL  +  Bn  cos  u>ntj 

n"-N /2 


T/N 


,  N  even 


(29) 


The  coefficients  An  and  Bn  are  computed  using  Equation  25. 
For  example,  let 


s 


G(s) 


a 

s  +  a 


so  that 


- 1  (GM)T/N 
N 


- 1  -  -  e~aT/N  (  1  - 

ff(z  .  e-aT/N)  (j  _  z-l] 


(30) 


(31) 


It  is  instructive  to  plot  the  Bode  plot  for  Equation  31  using  N  as  a 
parameter.  For  the  sake  of  clarity,  we  will  plot  versus  w  rather  than 
log  w  and  omit  the  phase  angle  plot.  Also,  for  reasons  of  clarity,  the 
ordinate  scales  will  be  displaced  for  the  different  values  of  N  (refer 
to  Figure  5).  Over  the  plotted  range  of  8ir,  the  N  *  1  case  repeats 
itself  4  times.  In  a  like  manner,  the  N  *  2  case  repeats  twice,  whereas 
N  =»  4  goes  through  one  cycle. 

Using  [1  sin  (n/2)t]  as  an  input,  in  the  N  »  1  case  our  only 
Interest  is  matching  the  sampling  points  with  a  single  sine  wave.  The 

magnitude  and  phase  angle  (not  shown  in  Figure  5)  could  be  read  from 

this  plot  at  w  -  tt/2,  it/2  +  2tt  ,  it /2  +  4tt  ,  it  /2  +  6n ,  •••;  giving  the 
correct  values  for  each  point.  Assume  next  that  the  input  has  a  fre¬ 
quency  b  =  (it/2  +  4tt).  Clearly,  if  the  objective  is  to  match  the  sample 

points  with  a  single  sinusoid,  the  frequency  of  the  output  could  be  b 

plus  any  2tt/T  multiple  as  the  sampler  cannot  tell  the  difference.  In 
fact,  the  "sub"  aliases  at  b  -  2w/T,  b  -  4n/T  will  also  work.  These 
"sub"  aliases  are  the  difference  terms  so  prominent  in  modulation 
theory. 

Figure  6  depicts  this  situation  for  a  steady-state  response  given  an 
input  frequency  of  b  ■  rr/2  rad/sec,  assuming  the  system  in  Equation  31. 


-  ----- 
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For  the  sake  of  clarity,  only  two  of  the  many  waves  which  fit  the  sample 
points  are  shown  —  one  at  w/2  rad/sec,  the  other  at  [(w/2)  -  (2n /T ) ] 
rad/sec • 

The  N  -  1  plot  in  Figure  5  corresponds  to  the  "sampled  spectrum"  of 
sampled  data  control  theory.  Turn  now  to  the  N  *»  2  case  wherein  the 
objective  is  to  match  one  inter-sample  point  as  well  as  the  sample 
points.  Let  the  input  frequency  be  it/2  and  note  that  the  points  at 
ui  =  ir/2,  tr/2  +  2n  give  the  correct  answers,  as  would  the  points 
n/2  +  4n,  n/2  +  6ir.  Suppose  next  that  the  input  frequency  is 
b  »  n/2  +  2n .  Clearly,  the  second  required  component  could  be  read  from 
the  "first  alias"  at  b  +  2n/T  or  the  first  sub-alias  at  b  -  2n/T  [or, 
for  that  matter,  all  the  frequencies  -  (n/2)  +  2mt  where  n  is  an  inte¬ 
ger]  . 

In  the  N  •  4  case,  four  sine  waves  are  required  to  fit  three  inter¬ 
sample  points  as  well  as  the  sampled  points.  If  the  input  frequency 

were  b  «  n/2  +  6»,  and  if  the  plot  of  Figure  5  with  its  limited  range 

of  8n  were  the  only  one  available,  clearly  it  would  be  to  our  advantage 

to  use  the  "difference"  frequency  points  «t  at  »  b  -  2*/T,  b  -  4n/T,  and 

b  -  6n/T  to  establish  the  magnitude  -d  relative  phase  of  the  three 

remaining  sine  waves. 

This  brief  discussion  serves  to  point  out  that  the  aliases  and  sub- 
aliases  can  be  associated  with  the  sum  and  difference  frequencies  of 
modulation  theory.  One  should  not,  however,  think  that  both  sum  and 
difference  components  must  be  simultaneously  present  in  the  output. 
Clearly,  only  N  components  are  needed.  We  can  now  remove  the  earlier 
constraint  that  b  <  2w/T.  If  b  is  less  than  2w/T,  it  is  certainly  true 
that 

“n  m  k  f  j  »  n  m  0,  1,  2,  ...»  N— 1  (32) 

However,  if  b  >  2w/T,  let  <a8  ■  2if/T  and  use 


>» 


to  restate  Equation  32  as 

%  -  b  +  ,  n  -  n0,  nQ  +  1,  ...»  0,  1,  2 . N  -  nQ  -  1  (34) 

For  our  previous  example  where  N  **  4,  suppose  b  ■  6it  +  ir / 2 .  Then 


(35) 


=  b  + 


2xn 

T 


-3,  -2,  -1,  0 


and  we  use  three  sub  aliases.  If  b  *  it / 2 ,  then 


n 


o 


2x  INX 


0 


(36) 


and  we  use 


wn  ‘  j  ♦  n  =  0,  1,  2,  3 
the  "positive"  aliases. 

Keep  in  mind  that  all  this  represents  a  convention  which  the  reader 

may  not  necessarily  elect  to  follow.  What  is  important  is  a  clear 

understanding  which  will  permit  one  to  pick  a  consistent  set  of  N  points 

from  the  Bode  plot.  In  this  regard,  the  reader  should  note  that  the  use 

of  Equation  28  (or  Equation  29)  instead  of  Equation  27  eliminates  the 

need  for  the  definition  of  n  . 

o 

Of  interest  is  the  case  where  N  is  extremely  large.  In  fact,  let 
N  *  08  after  evaluating  Equation  31  at  z  ■  1  iui^T/N): 


21 


z«14%(T/N) 


11m 


(37) 


■  -X.^-aJ/N  ... 

N(z  -  e-aT^N )  1  -  z_1 


(l  -  e_aT/N)(l  -  U-unT) 


NlUWnfT/N)  -  e-*T^][l  -  U-0)n(T/N)] 


An  indeterminate  form  is  obtained.  Therefore,  use  L'Hopital's  rule 
twice  (substitute  liui^T/N)  •  cos  (^(T/N)  +  J  sin  (^(T/N),  etc.)  and 
obtain: 


/l  *-s_T\|T/N 

-In)  T 

1  -  e  J  n  1 

1  -  e~sT  1 

U(s  +  l)/|z-l*wn(T/N)  |11(n 

j“nT  1  +  J“n 

sT  s  +  1 

n-» 


(38) 


That  is,  as  N  ■*  one  simply  divides  GM  by  T  and  evaluates  the  coeffi¬ 
cients  at  s  »  ju>n.  This  is  representative  of  the  general  result  dis¬ 
cussed  in  the  next  subsection.  A  word  of  caution  is  in  order  on  <un.  As 
N  approaches  infinity,  the  definition  of  Equation  28  (or  Equation  29) 
should  be  adhered  to  in  order  to  properly  incorporate  the  sub-aliases 
into  the  spectrum. 

F.  OPEN  LOOP  FREQUENCY  RESPONSE  — 

CONTINUOUS  OUTPUT 

The  next  development  uses  an  approach  suggested  by  Professor  Hsi-Han 
Yeh  (University  of  Kentucky) .  This  approach  more  clearly  shows  the 
dependence  of  the  spectrum  on  both  positive  and  negative  aliases.  That 
is,  for  finite  N,  there  is  a  choice  in  the  makeup  of  the  sinusoidal  com¬ 
ponents  which  exactly  match  the  steady-state  sample  points.  The  finite 
N  case  is  therefore  not  unique  —  in  sharp  contrast  to  the  infinite  N 
results  which  require  (as  will  be  shown)  the  use  of  all  the  sub-aliases 
as  well  as  all  of  the  aliased  components. 
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In  the  previous  section  it  was  shown  that 


z*esT/N 

An  +  JBn  *  n  (GM>T/N|  (39) 

z*Uun(T/N) 


To  deduce  the  behavior  for  infinite  N,  use  an  alternative  description  of 
Figure  4  which  describes  the  sine  wave  input  directly  in  the  s-domain: 


C(s) 


G(s)M(s) 


_1 

T 


s^  +  wr 


(40) 


where  =  b  +  2irn/T. 


The  partial  fraction  expansion  of  Equation  40  may  be  written  as 


C(s) 


Ant^n)  Bn3 


+  [Transient  modes  of  GM] 


(41) 


=-“>  S^+dlj!  Sz+0)n 

Multiply  each  side  of  Equation  41  by  s2  +  and  evaluate  at  s 
C  ( s )  ( s 2  +  uj)  =  +  Bnja)n 


3V 


-  G(s)M(s)  (1  /T)u>n 


An  +  JBn  *  (l/T)G(s)M(s) I  ,  n  =  0,  +1,  +2,  ...  (42) 

n  n  |s=jwn 

Thus  the  continuous  spectrum  contains,  because  of  the  summation 
from  -*•  to  +°°,  both  positive  and  negative  aliased  frequencies. 

The  finite  N  example  of  the  previous  section  can  now  be  studied  for 
the  case  of  infinite  N.  Thus 


An  +  IK  “ 


l---.a-3T  . 

sT 


+  1 


s-ju>n 


produces  the  Bode  plot  of  Figure  7. 


(43) 
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The  interpretation  of  Figure  7  is  as  follows.  Suppose  a  unit  sine 
wave  at  1  rad/sec  is  input  to  the  sampler.  Then,  if  sine  waves  at  1, 
1  +  2n/T,  1  -  2x/T,  1  +  4x/T,  1  -  4*/T,  ...  are  added  together,  the 
resultant  waveform  will  be  an  exact  match  of  the  actual  steady-state 
output  waveform.  In  Figure  7,  one  may  plot  the  sub  aliases  (the  nega¬ 
tive  frequency  components)  on  a  "positive  frequency"  Bode  plot  by  taking 
advantage  of  the  fact  that  the  magnitude  is  an  even  function  of  fre¬ 
quency  and  the  phase  is  an  odd  function  of  frequency. 

One  would  expect  this  waveform  to  be  relatively  clean,  since  the 
first  alias  is  30  dB  lower  than  the  input  component.  However,  the  tran¬ 
sient  response  itself  does  not  bear  out  this  conjecture,  as  can  be  seen 
in  Figure  8.  The  reason  is  that  the  higher  terms  are  important.  They 
do  not  represent  "harmonic"  terms  but  are  rather  modulation  components 
which  must  add  together  properly  in  order  to  match  conditions  at  the  T 
transition  points.  It  can  be  seen  that  the  "steady  slate"  does  not 
necessarily  take  on  the  additional  attribute  of  periodicity.  This 
occurs  only  when  the  input  frequency  and  the  sampling  frequency  bear  an 
integer  relationship  with  respect  to  one  another. 

G.  INPUT  SIGNAL  WITH  PHASE  SHIFT 

If  the  Input  signal  has  the  form 


r(t)  -  kj  sin  bt  +  k2  cos  bt  (44) 

the  results  of  the  previous  section  are  changed  only  by  a  complex  con¬ 
stant.  Following  exactly  the  procedures  of  Subsection  E,  except  for 
using  the  more  general  input  given  in  Equation  44,  gives  (Reference  1): 


2£e8T/N 

An  +  jBn  -  (GM)T/n|  .  (kl  +  jk2) 

z-Hwn(T/N) 


(45) 


Given  that  the  limit  N  +  «: 


An  +  JBn  «  4  CM 


s-ju>n 


•  (kj  +  jk2) 


(46) 
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Figure  8.  "Steady-State"  Transient  Response 

The  choice  of  the  N  frequencies  for  the  finite  N  case  is  at  the  discre¬ 
tion  of  the  user.  However,  both  "negative"  and  "positve"  frequencies 
are  required  as  N  +  »  (Equations  28  and  29). 

H.  SINGLE-RATE  CLOSED-LOOP  FREQUENCY  RESPONSE 

The  closed-loop  results  will  be  configuration  dependent.  However, 
the  mathematics  remains  tractable  and  can  be  followed  through  on  a  case- 
by-case  basis.  It  is  important  to  have  an  insight  into  the  mathematical 
structure  and  the  particular  simplifications  that  surface  in  a  closed- 
loop  analysis. 

Consider  the  (vector)  system  shown  in  Figure  9.  The  procedure  we 
now  follow  will  be  typical.  First,  solve  for  the  vector  component  at 
the  input  of  the  data  holds. 

Et  -  g]rT  -  GjG 2 (GM) TET  (47) 

Therefore 

Et  -  [l  +  g]g2(GM)t]  GiRt  (48) 
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Figure  9.  Illustrative  Vector  Closed-Loop  Configuration 
Next,  solve  for  C(s): 

C  =  ( GM ) [ I  +  GiG2(GM)t]  ^gJrt  (49) 

The  spectrum  of  C(s)  is  of  interest;  we  seek  it  by  finding  first  the 
spectrum  of  C^^  and  going  to  the  limiting  case  of  N  +  ». 

Let  the  input  be  a  sine  wave  at  frequency  b  rad/sec  and  let  the 
delay  operator  be 

Z-1  =  e-sT/N  (50) 

so  that 

RT 

Therefore 

CT/N  *  [(GM)[I  +  GiG2(GM)t]-1g]rT]T/,N 

CT/N  =  (GM)T/N[!  +  GiG2(GM)T]  ^[r1 


For  the  sake  of  brevity  write  Equation  52  as 


sin  bT 


Z2N  ^  2  [cos  bT )z^  +  1 


(51) 


Note  that  Equation  53  is  exactly  the  same  as  Equation  13,  except 
T /N  T /N  T 

(GM)  has  been  replaced  by  Gg.  Hence  we  can  write  a  key  result 
using  Equation  25: 


But 


An  +  JBn 


1  T/N  T 

N  "  B  z-H<un(T/N) 


Gj  ±  Gb(zn) 


(54) 


Therefore,  using 

[H<VT/N)]n  “ 

*  cos  wnT  +  j  sin  u)nT 

-  cos  [b  +  (2WT)]T 

+  j  sin  [b  +  (2nn/T) ] T 

“  cos  bT  +  j  sin  bT  (55) 

we  obtain 


Gb  14^(T/N) 


G[Uwn(T/N)]N 
G[UconT]  -  G  (1  4  bT] 


(56) 


N  T 

It  is  permissible  to  replace  z  in  Gg  with  z  and  evaluate  it  at 
z  ■  l4bT.  At  this  point  we  have 


An  +  JBn 


a  sT/N 
z*e 

z-lAu^CT/N) 


Gg(z) 


,  si  \ 

z-e  \ 

z*»  l4bT/ 


(57) 


Equation  57  is  the  basic  result  for  the  finite  N  case.  To  reiterate,  to 
find  the  coefficients  of  the  N  sine  waves  for  the  T/N  sampled  output  of 
C,  compute  the  normal  "T"  transfer  functions  for 


T  T  T  —1  T 

Gg  -  [I  +  G1G2(GM)TJ  Gj 


(58) 


and  evaluate  It  for 


z  -  UbT 


(59) 


Next,  compute  the  normal  T/N  pulsed  transfer  function  for  and  evalu¬ 
ate  it  at  z  *  l4Wn(T/N)  where  un  “  B  +  (2itn/T) . 

T 

Thus,  Gu  is  periodic  in  (2fln/T)  and  it  suffices  to  use  bT  instead 

T/N 

of  u>n(T/N) .  Moreover,  only  the  is  a  function  of  N;  this  simplifies 

the  procedure  involved  in  the  limiting  case  tremendously.  For  the  case 
of  N  ♦  the  continuous  case,  we  obtain: 


^  +  JBn  = 


M(s)G(s) 


Gj(z) 


C  cjT 
z-e13 1 


z=  I4  bT ) 


0,  +1,  +2, 


(60) 


Equation  60  is  the  desired  result  for  the  given  closed-loop  configu¬ 
ration.*  However,  the  mathematical  ideas  are  what  count;  one  can  follow 
the  details  through  for  other  configurations  with  relative  ease. 

With  this  development,  one  is  in  a  position  to  plot  the  Bode  plots 
for  the  closed-loop  system  of  Figure  1  and  verify  the  results  given  in 
Fig.  2.  We  will  also  use  that  example  to  solidify  the  meaning  of  the 
frequency  response  for  the  finite  N  case.  Suppose  the  continuous 
transient  response  for  a  unit  amplitude  sine  wave,  with  a  frequency 
of  tt/2,  is  available  (see  Fig.  10).  According  to  theory,  one  should  be 
able  to  set  N  =  l  and  from  the  Bode  plot  read  the  magnitude  and  phase  of 
the  single  sinusoid  that  fits  the  sample  points  at  the  sample  instant. 
This  indeed  proves  to  be  the  case  and  is  shown  in  Figure  11.  In  Fig¬ 
ure  11  a  section  of  the  transient  response  has  been  "copied"  and 


Accurate  numerical  determination  of  Gglj^bT  nay  prove  difficult  at 
high  sampling  rates.  This  is  the  result  of  small  differences  between 
large  numbers  which  occur  in  the  computations  as  poles  and  zeros 
approach  the  unit  circle.  In  this  event,  one  is  well  advised  to  carry 
out  equivalent  computations  in  a  domain  where  numerical  conditioning  is 
much  improved  (e.g.,  in  terms  of  w'  or  w) . 
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Steady-State  Sinusoidal  Components 


overlaid  with  the  sine  wave  which  results  from  the  N  ■  1  computation 
namely 


(C(t)]T 


7T 

sin  — 


,  71  1 
L  +  b0  cos  —  t  J 


(61) 


Equation  61  indicates  that  the  value  of  the  wave  at  the  sampling 
instants  is  of  interest.  However,  for  expository  purposes,  a  complete 
cycle  has  been  shown. 

Next,  consider  the  N  =  2  case  where  the  desire  is  to  fit  not  only 
the  sample  point  but  one  inter-sample  point  as  well.  This  case  is  shown 
in  Figure  12.  The  T/2  response  equation  is 


C(t) 


T/2 


ss 


IT  TT 

aQ  sin  —  +  bg  cos  —  t 


+  a>  sin 


TT  2  IT  "I 
9  +  “H  +  bl 


(2 L  .  2TTV 
cos  [-  +  — Jt 


T/2 


(62) 


Again,  for  expository  reasons,  the  continuous  waveform  is  shown  which 
results  from  sines  and  cosines  at  cu  =  n/2  and  its  first  alias  at 
w  *  (tt/2)  +  (2tt/T).  A  half  period  for  the  N  <=  10  case  is  shown  in  Fig¬ 
ure  13.  The  steady-state  wave  of  this  example  is  periodic  and  free  of 
modulation  effects  simply  because  the  selected  input  frequency  bears  an 
integer  relationship  to  the  sampling  frequency. 

It  is  also  important  to  bear  in  mind  when  "matching"  sample  points 
that  the  Bq  term  for  N  =  1  will  be  unequal  to  the  Sq  term  for  N  ■  2. 
This  is  demonstrated  in  Table  l  for  N  =  1,  2,  and  4  and  b  *=  n/2. 


I,  A  PARTICULAR  TWO-RATE  CONFIGURATION 


As  with  the  analysis  of  closed-loop  single-rate  systems,  the  analy¬ 
sis  of  the  multi-rate  closed-loop  case  is  configuration  dependent.  Here 
we  treat  a  particularly  simple  two-rate  configuration  which  can  be  ana¬ 
lyzed  without  resorting  to  a  switch  decomposition  format.  A  more  com¬ 
plex  example,  which  requires  the  use  of  vector  switch  decomposition, 
will  be  treated  after  the  fundamentals  of  that  technique  have  n 
reviewed  in  Section  III. 


TABLE  1.  COMPONENT  COEFFICIENTS  N  =  1,  2,  A 


N  -  1 

N  =  2 

N  =  A 

b 

aQ  -0.17AA68021 

bQ  -0.2876A9137 

aQ  -0.0A8579760 

b0  -0.306381976 

aQ  0.002998683 

b0  -0.302606086 

b 

a:  -0.125888261 

bL  0.018732839 

a2  -0.0A6197029 

b2  -0.0A35A8018 

a2  -0.051578AA3 

b2  -0.003775889 

a3  -0.079691232 

b3  0.062280857 

Consider  the  two-rate  system  shown  in  Figure  1A.  In  Figure  1A,  Wj 
and  are  compensation  networks,  M  is  a  data  hold,  and  G  represents  the 
open-loop  system  dynamics  (consider  these  to  be  matrices  of  the  proper 
dimensions ) . 

The  objective  is  to  find  the  "frequency  response"  for  the  output 
vector  C.  As  in  the  single-rate  case,  assume  that  C  undergoes  a  phantom 
T/N  sampling  operation  and  then  seek  the  limit  as  N  +  »,  From  Fig¬ 
ure  1A: 

C  =  GMET/M  (63) 

or 

CT/N  „  (gm)T/NET/M  t  N/M  an  integer  (bA) 

The  first  task  is  to  solve  for  E^^M.  This  is  non-trlvial;  the  details 
must  be  followed  with  care. 

E  -  W iRt  -  WiW^GMET/M)1  (63) 
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Figure  14.  A  Specific  Two-Rate  Closed-Loop  Configuration 
Therefore, 

eT/M  «  -  W 1^2  (GMET/m)  (66) 

Pre-multiply  Equation  66  by  GM  and  sample  at  a  T  interval. 

(GMET/M]T  *  (GMwY/M3  RT  -  (GMwJ/M)  W2(Q1Et/m)T  (67) 

Solve  Equation  67  for  (GMET/M)T: 

(GMET/M)T  -  [i  +  (GMWi/M)TW2]  (gMw|/M)TRt  (68) 

Substitute  Equation  68  into  Equation  66  and  clear  through.  The  result 
is : 

eT/M  „  WT/M{ x  _  y\[i  +  (gMWi/M)  W*]  (GMWi/M)  }rt  (69) 


For  brevity,  let 


et/m  -  Wi/mgJrT 


T  (  T  /  H  ^ 

The  evaluation  of  GA  ia  not  elementary.  For  example,  the  [GMWi  )  ele- 

T 

ment  of  GA  will  have  to  be  computed  using  either  switch  decomposition  or 


the  phantom  sampler  (Reference  1  or  see  the  development  of  Section  IV) 
Via  the  phantom  sampler, 


(GMWi/M)  =  [(GM)T/Mwj/M] 


(71) 


To  this  point,  the  two-rate  example  yields 
CT/N  _  (q<)T/NwJ/MgIrT 


(72) 


and  we  see  that  the  only  new  element  added  over  the  single-rate  case  Is 
the  addition  of  a  term  sampled  on  a  T/M  Interval  together  with  a  con¬ 
straint  that  N/M  be  an  integer. 


Let 


z  -  e 


sT/N 


(73) 


so  that  a  unit  amplitude  sinusoidal  input  at  b  rad/sec  has  the  transform 


RT 


_ S-N.  AiSL  fr.I _ 

Z2N  _  (2  cos  bT)zN  +  1 


(74) 


As  in  the  single-rate  case. 


substitute  Equation  74  into^Equation  72: 


CT/N 


(GM)t/nWj/MgJ 


_ zN  sin  bT _ 

z2N  -  [2  cos  bl)zN  +  1 


(75) 


Again,  the  problem  is  in  a  recognizable  form  and  we  may  proceed  directly 
to  the  answer. 


*»  +  ’  »  (GM)T/NW^M<:*jl,.„(T/N)  (76> 

Next,  since  the  local  definition  of  z  ■  is  in  effect,  express 

e®T  .  [esT/Nj^  and  vrite 


Therefore 
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ga<*n>  |  u^n/N)  - 

-  GA(l4«)nT)  -  GA(l4bT) 


(77) 


That  is,  take  the  "T"  2-transform  of  GA,  and  evaluate  at  z  •  hJ-bT.  Now, 
the  "new"  element, 

W^/M  A  w^zN/M)  (78) 

Therefore, 

“1<zN/M||i,«„(T/H)  ■  »l[U«n(T/N)]N/M  m) 

-  W1[l5.0)n(T/M)]  (80) 

taking  the  "T/M"  z-transform  of  Wj  and  evaluating  it  at  z  -  Kwn (T/M) . 

At  this  point,  only  GM  depends  on  N  and  we  can  go  to  the  limit  of 

N  -*■  «. 


An  +  J®n 


J. 

N 


ziesT/N 


(GM)T/N 


z-UUn(T/N) 


wT/M 

W1 


zlesT/M 

z«esT 

_T 

ga 

z-lAu>n(T/M) 

z-l$bT 

(81) 


This  is  the  desired  result  for  finite  N  (remember  N  -  M,  2M,  etc.)* 
As  N  >  the  coefficients  of  the  continuous  spectrum  are  given  by 

A„  +  jBn  -  ^  (GM) 

1  a-J<% 

n  ■  0,  +1,  +2,  ... 


2AeaT/M 

-  sT 
z«e 

wf" 

gt 

gA 

z-l^u^T/M) 

z-UbT 

(82) 
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The  algebraic  manipulations  used  in  this  example  circumvented  the  need 
for  the  use  of  a  switch  decomposition  model  by  using  the  classic  concept 
of  a  phantom  sampler*  The  limitations  of  this  type  of  algebraic  manipu¬ 
lation  will  be  discussed  in  more  detail  in  Section  V. 

J.  SECTION  SUMMARY 

The  "sampled  spectrum"  concept  of  sampled  data  control  theory  is 
concerned  with  determining  the  single  sinusoid  which  fits  the  output 
samples  of  a  single-rate  system  at  the  sampling  instants.  In  this  sec¬ 
tion,  we  have  reviewed  the  extension  which  encompasses  the  continuous 
spectrum  of  the  continuous  variables  in  a  discretely  controlled  system. 
Moreover,  the  theory  considers  the  finite  N  case  wherein  one  is  con¬ 
cerned  with  the  group  of  N  sinusoids  that  matches  the  data  not  only  at 
the  sample  points  but  at  N-l  inter-sample  points  as  well.  This  is  an 
important  aspect  since  bench  validation  of  digital  hardware  is  often 
specified  in  terms  of  an  end-to-end  "frequency  response."  Since  output 
data  are  taken  at  a  finite  number  of  points,  it  will  be  important  to 
compute  finite  N  results;  the  coefficients  may  differ  significantly  from 
the  continuous  (N  +  00 )  values. 

The  results  for  the  closed-loop  cases  have  been  configuration  depen¬ 
dent;  however,  the  basic  technique  is  relatively  clear.  One  starts  at 
the  continuous  state  vector  and  writes  the  system  equations  back  to  the 
first  input  point.  The  next  fundamental  step  is  to  convert  that  input 
into  an  equation  which  contains  the  sinusoidal  input  as  the  basic  forc¬ 
ing  function.  At  this  point,  the  basic  equations  which  apply  to  the 
open-loop  case  can  be  used.  A  closed-loop  two-rate  case  was  discussed 
which  did  not  require  the  use  of  switch  decomposition.  The  limitations 
of  the  algebraic  manipulations  will  be  discussed  more  fully  in  Sec¬ 
tion  V. 
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SECTION  III 


VECTOR  SWITCH  DECOMPOSITION  AND 
A  "SCALAR"  APPROACH 


A.  INTRODUCTION 

This  section  reviews  vector  switch  decomposition,  noting  the  dimen¬ 
sionality  problems  introduced  by  multiple  frame  times.  This  is  followed 
by  a  brief  discussion  pertaining  to  a  class  of  multi-rate  problems  which 
can  be  treated  without  recourse  to  switch  decomposition.  Several  exam¬ 
ples  are  then  used  to  demonstrate  the  scalar  algebraic  manipulation 
required,  establish  the  relationship  between  the  scalar  algebra  and 
switch  decomposition,  and  motivate  the  need  for  an  algorithm  which  can 
be  used  to  evaluate  the  various  expressions  that  result.  The  algorithm 
itself  is  discussed  in  Section  IV. 

B.  REVIEW  OF  SWITCH  DECOMPOSITION 

In  essence,  switch  decomposition  is  a  procedure  wherein  systems  hav¬ 
ing  multiple  sampling  operations  (occurring  at  fixed  but  unequal  sampl¬ 
ing  intervals  but  with  a  sampling  pattern  which  is  repeated  over  a 
fixed,  finite  time  interval)  are  converted  into  an  equivalent  single 
sample  rate  format.  As  originally  Introduced  by  Kranc  (Reference  3), 
the  method  used  a  summing  point  methodology  that  proved  to  be  extremely 
cumbersome  when  the  ratios  of  the  sampling  periods  become  high.  For 
this  reason,  and  because  evolving  state  transition  methods  were  tending 
to  push  transform  methods  into  the  background,  the  method  fell  Into  dis¬ 
use.  However,  there  is  much  to  recommend  the  switch  decomposition  con¬ 
cept  for  use  in  both  time  domain  and  transform  domain  analyses.  In  the 
subsection  that  follows  we  will  review  the  basic  concept  and  remove  some 
earlier  restrictions  by  recasting  It  in  vector  form.  The  vector  form 
simplifies  matrix  block  diagram  manipulation  for  multiloop,  multi-rate 
sampled  systems. 


An  example  reveals  the  key  ideas.  Consider  the  continuous  signal, 
shown  in  Figure  15a,  to  be  sampled  at  T/3  samples  per  second.  This 
results  in  the  sample  sequence  shown  in  Figure  15b.  The  sampled  values 
have  been  numbered  for  easy  reference.  Suppose  we  now  sample  the  con¬ 
tinuous  signal  with  a  sampling  period  T.  This  results  in  the  sample 

sequence  consisting  of  sample  numbers  1,  4,  7,  10,  13,  ...  shown  in 

T 

Figure  15c.  Define  this  sample  sequence  to  be  R  . 

Next,  advance  the  continuous  signal  R  by  T/3.  Then  sample  the 
advanced  signal  with  a  sampling  period  T.  This  results  in  a  sample 
sequence  consisting  of  sample  numbers  2,  5,  8,11,  14,  ...  shown  in  Fig¬ 
ure  15d.  Define  this  sample  sequence  to  be  (esT^R)T.  Finally,  advance 
the  continuous  signal  R  by  2T/3  and  sample  it  with  a  sampling  period  T. 
This  results  in  the  sequence  consisting  of  sample  numbers  3,  6,  9,  12, 
15,  ...  shown  in  Figure  15e.  Define  this  signal  sequence  to  be 

(=2sT/3R)T. 

The  significance  of  the  switch  decomposition  concept  resides  in  its 

ability  to  provide  an  alternative  expression  for  the  original  sequence 
T/3 

R  in  terms  of  several  quantities  which  are  each  sampled  simultane- 

T/3 

ously  every  T  seconds.  This  alternative  expression  for  R  consists  of 
the  sum  of  RT,  (esT^R)T,  and  (e^sT/,3R)T: 

rT/3  ,  RT  +  (esT/3  r]T  e-sT/3  +  (e2sT/3  R)T  e-2sT/3  (83) 

Equation  83  has  a  simple  factored  equivalent  that  is  the  product  of  two 
vectors  and  the  scalar  R, 


_T/3  ,  -sT/3  -2sT/3 

R  -  1,  e  ,  e 


Or,  more  compactly, 


,T/3  „ 


where 


and 


W 


W(W*R); 


W  =  [l, 


-sT/3  -2sT/3l 


1 

esT/3 

2sT/3 


C.  EXTENSION  TO  THE  VECTOR  CASE 


(85) 

(86) 

(87) 


Further  generalization  allows  R  to  be  a  vector  of  continuous  sig¬ 
nals.  It  is  necessary  to  define  a  least  common  sampling  period,  T,  and 
a  greatest  common  subinterval,  T0»  respect  to  the  R  vector.  The  p 

elements  of  R  may  be  sampled  at  different  minor  sampling  periods:  T^, 

. . . T  ,  respectively.  It  is  further  assumed  that  the  minor 

sampling  periods  are  such  that  a  finite  positive  T  exists  such  that 


T  -  NjTi  -  NiTi  •”  =  NpTp  (88) 

holds  for  a  set  of  finite  positive  integers: 


N i ,  • • * »  ,  •  •  • ,  Np 

The  minimum  T  for  which  Equation  88  holds  is  the  least  common  sampling 
period  (for  R) .  A  subinterval  can  be  found  for  which 

T  -  NT0  (89) 

and  N/N^  is  an  integer  for  all  i  *  1,  2,  ...,  p.  The  largest  value 

of  TQ  satisfying  these  conditions  is  the  greatest  common  subinterval 
(for  R) .  Equation  89  defines  N  for  the  greatest  common  subinterval. 
Given  values  for  N,  N^,  p>  and  T,  the  p  *  block  diagonal  matrix,  W, 


0 


W  -  W(s) 


where 


(90) 


Wl  _  [i  e-sT/Ni  ...  e-s(j-l)T/Nl  ...  ^sCNi-DT/Nij 


(91) 


The  operator  matrices  W  and  W  can  be  used  to  represent  multi-rate 
sampling  operations  in  terms  of  a  single-rate  sampling  operator  in  vec¬ 
tor  block  diagrams.  This  is  illustrated  in  Figure  16. 

Consider  an  example.  Let  R  be  a  vector  with  Lhree  components.  Let 
the  first  component  be  sampled  with  period  T/6,  the  second  with  period 
T/3,  and  the  third  with  period  T/2,  i.e., 


R 


* 


(92) 


The  objective  is  to  compute  W  in  order  to  obtain  an  explicit  expression 

* 

for  R  via  Equations  90  and  91  (which  is  equivalent  to  Figure  16).  For 
this  example. 


p  *  3 

T  is  the  least  common  sampling  period 
Tj  -  T/6,  T2  -  T/3,  T3  -  T/2 

N1  “  6>  n2  “  3»  N3  "  2 

T/6  is  the  greatest  common  subinterval 

N  -  6 
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Therefore 


1,  e-sT/6,  e-2sT/6f  e-3sT/6f 

e-4sT/6>  e-5sT/6 

0 

0 

w  - 

i 

I 

o 

e-2sT/6>  e-4sT/6 

0 

1 

0 

0 

lf  e-3sT/6 

(93) 

This  example  gives  some  insight  as  to  how  increased  dimensionality  can 

complicate  problems  in  practical  application.  Consider  the  vectors  R, 
*  T 

R  ,  and  (W*R)  .  These  vectors  have  3,  3,  and  11  elements,  respectively. 
The  vector  (w*R)T  will  have  pEN^  elements  in  general;  whereas  R  and  R* 

will  have  only  p  elements  each.  This  is  significant  in  that  analyses 

T 

will  tend  to  be  conducted  in  terms  of  vectors  like  (W*R)  in  distinction 

* 

to  vectors  like  R  •  Consequently,  the  potential  for  expanded  dimension- 
ality  in  connection  with  analyses  of  multi-rate  sampled  problems  is 
great.  For  example,  consider  a  problem  wherein  there  are  two  minor 
sampling  periods,  39  ms  and  40  ms.  It  is  easy  to  verify  that  the  dimen¬ 
sionality  expansion  factor,  N,  is  1560. 

On  the  more  positive  side,  matrix  operations  are  routine.  Consider 
the  system  shown  in  Figure  17.  Once  the  vector  multi-rate  sampling 
operations  in  Figure  17a  have  been  replaced  by  the  switch  decomposition 
equivalent  (Figure  17b),  analytical  manipulations  are  routine: 

E1  “  (WUR)T  "  (WUHW2)T(W2aGW1)T  Ei  (94) 
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D.  NONSYNCHRONOUS  SAMPLING 


Nonsynchronous  sampling  is  a  basic  tool  useful  for  modeling  distrib¬ 
uted  computation  architecture,  data  skevmess  in  the  A/D  and  D/A  conver¬ 
sion  processes  as  well  as  the  internal  computational  delay  of  the  digi¬ 
tal  computer-  By  definition,  nonsynchronous  sampling  occurs  when  all 
the  systems'  sampling  operations  are  repeated  at  the  same  rate  but  occur 
at  different  instants  of  time  (refer  to  Figure  18). 

In  Figure  18  both  continuous  signals,  Xj  and  X2»  are  sampled  at 

1/T  Hz,  but  the  X2  sampler  is  "out  of  sync"  with  the  x^  sampler  by 

T 

TQ  seconds.  The  sampling  operation  for  x*  is  shown  symbolically  in  Fig¬ 
ure  19a  and  for  x2  in  Figure  19b.  (*  notation  on  X2  is  used  here  to 
indicate  an  "unconventional"  sampling  operation.) 

Figure  19b  models  the  nonsynchronous  sampler  with  a  synchronous 
sampler  by  preceding  the  sampler  with  the  operation  followed  by  the 
operator  W,  i.e., 

-  W(W*W2)T  (100) 

where 

-sT  sT 

W-e°  ,  W*«e°  (101) 


Proceeding  according  to  Equation  101,  one  advances  X2  by  T  seconds, 
samples  at  the  1/T  rate,  and  then  delays  (WAx2)  by  TQ  seconds  to  obtain 
the  time  sequence  (refer  to  Figure  20).  Note  how  the  nonsynchronous 
sampling  operation  on  X£  is  modeled  in  terms  of  a  scalar  factor;  thus 
the  dimension  of  the  equivalent  single-rate  sampled  signal,  (W  X2)T,  is 
not  increased. 


The  model  readily  extends  to  the  case  where  x  is  a  vector 


E.  SPECIAL  CASE 


Consider  Figure  21,  a  three-rate  system  without  feedback. 

CX/»3  .  (g2(GiR^»1)I/N2JT/1‘3 

-  0^3  Gj/»2  RT/"i 


(102) 


(103) 


if  and  are  integers.  However,  Equation  103  can,  using  the 
algorithm  given  in  Section  IV,  be  evaluated  directly  (for  a  given  R)  in 
the  event  that  anc*  N2^Nl  do  not  sat-isfy  ihe  integer  relationship. 
This  is  not  the  case  for  the  feedback  system  of  Figure  22. 


E 


(104) 


If  M/N  is  an  integer,  Equation  104  becomes 


E 


R 


Gjof"  ET'" 


which  indicates 


ET/N  .  RT/N  _  [g2gT/M]T/N  ET/N  (105) 

Solve  for  the  unknown 

eT/m  .  [1  +  (g2gT/h)t/n]_1  rt/n 

using,  for  example,  the  algorithm  given  in  the  next  section.  If  M/N  is 
not  an  integer,  one  would  be  forced  to  cast  Figure  22  in  a  vector  switch 
(or  some  other  equivalent  state  transition)  format. 


i 


Figure  22.  A  Closed-Loop  Configuration 


In  a  like  manner,  if  N/M  is  an  integer,  write 
c  -  gJr  -  G2Ct/m]T/N 
-  glrT/n  .  GigT/nct/m 

giving 

CT/M  -  [glrT/N]T/M  -  [giGt/n]T/m  ct/m 

GT/M  »  [i  +  (g^G^/^)^^M]  [GjR^/n]^M  ,  N/M  is  an  integer 


(106) 


To  summarize,  open-loop  systems  such  as  Figure  21  can  in  theory  be 
analyzed  without  recourse  to  switch  decomposition.  Closed-loop  systems. 
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on  the  other  hand,  require  a  fortuitous  set  of  relationships  among  the 
frame  times. 

F.  A  PARTICULAR  THREE-RATE  CLOSED-LOOP  SYSTEM 

There  is  a  natural  tendency  in  the  design  of  digital  systems  to 
select  sampling  ratios  in  powers  of  2.  This  is  probably  due  to  the  ease 
with  which  clock  rates  can  be  doubled  (or  halved).  Such  is  the  case  of 
the  three-rate  digital  controllers  for  the  Shuttle  (25,  50,  100  Hz)  and 
the  F-18  (20,  40,  80  Hz).  It  has  been  our  experience  that  such  systems 
can  always  be  analyzed  with  operations  such  as  those  described  in  the 
previous  section.  For  example,  consider  the  closed-loop  configuration 
of  Figure  23.  One  may  verify  the  following  operations: 


E  =  R  -  G2[G1Et/2]T  -  G3[G1ET/2]T/4  (107) 

E  =  R  -  G2[G1Et/2]T  -  G3G?/4  eT/2 

eT/2  =  RT/2  _  GT/2[GieT/2]T  -  [G3GTA]T/2  E?/2  (108) 


Therefore, 

ET/2  =  [i  +  (G3GT/4)T/2  j-1  1rT/2  _  GT/2(GiET/2)Tj 

-  [  •  ]_1  Rt/2  -  [  •  J-1  G|/2[G1ET/2]T  (109) 

where 

l-l  •  U  +  (C3GT/4)T/2] 

Equation  109  defines  ET^2  in  terms  of  R^/2  and  [GiE^/2]t.  Next,  multi¬ 
ply  Equation  109  by  Gj  and  solve  for  [GiE^/2]^: 


[GlET/2jT  -  [Gl(  •  ]-l  rT/2]T.  [Gl(  .  j-1  GT/2]T[GlET/2jT 


G 


G 


Figure  23.  A  Three-Rate  Configuration;  T,  T/2,  T/4 


[GiET/2]1  -  {1  +  [G1[*r1G2/2]TJ  [g1  [• )  ~1rT/2  ]T 


(110) 


Substituting  Equation  110  into  Equation  109,  and  using 

C  -  GjE^/2 

gives  the  continuous  output  in  terms  of  R.  One  may  contrast  the  com¬ 
plexity  of  this  algebra  against  the  simplicity  of  switch  decomposition. 
Since 


C  -  G^Ei  -  GxE1/2  -  Et/2 


dll) 


one  need  solve  only  for  e|: 


[I  +  (W2*G2)T(G1W2)T  +  (w2*G3W4)T(w4jtG1W2)T]El  -  (w2*r)T 


(112) 


Thus,  a  simple,  straightforward  block  diagram  manipulation  gives  the 
T 

needed  vector,  Ej,  essentially  "by  inspection": 

El  "  [I  +  +  (w2jkG3W4)T(w4jkG1W2)T]_1(w2jtR)T  (113) 


'■Hi'-? 


The  switch  decomposition  result.  Equation  113,  can  be  used  directly 
to  arrive  at  Equation  108.  Premultiply  Equation  112  by  and  "nest" 
the  four  terms*  Two  terms  are  obvious: 

W2eT  &  eT/2 

W2(W2*R)T  *  Rt/2 

The  other  two  are  slightly  more  laborious.  For  example: 

T  ( I  T|T 

(w2*G3w4)T(w4*GiW2)TE^  -  j[(w2AG3W4)(w4jkG1W2)E^J  j 

-  J[(w2*G3W4)(w4jkG1ET/2)]TJ 

(  r  t/21T/^T 

-  \v2tp3[Gl*V2]  j 

-  [m2*G3g{/4eT/2] 

but 

w2[w2*G3gT/AEt/2]  i  [G3Gi/4ET/2]  /2  .  [g3gJ/4]  '  ET/2  (H4) 

A  similar  operation  on  W2(w2aG2]t(GiW2)tEi  gives  the  remaining  terms 
of  Equation  108.  Specifically, 

W2(w2*G2JT(GiW2)te]'  =  Gr2tGlET/2]T  (115) 

It  can  be  appreciated  that  switch  decomposition,  coupled  with  the 

nesting  operation,  can  be  used  (where  appropriate)  to  generate  a 

"scalar"  alternative  to  the  switch  decomposition  model.  The  use  of  the 

T 

word  scalar  is  in  the  sense  that  E^  is  a  2  x  1  vector  in  Equation  113 

T/2  r  T/4i^'2 

whereas  E  in  Equation  108  and  [G3Gi  J  in  Equation  110  are  sca¬ 
lars.  The  use  of  the  "scalar"  results  implies  the  existence  of  an 

m  M 

algorithm  for  evaluating  terms  like  [g.jG!'*]  .  Just  such  an  algorithm 

la  given  in  Section  IV. 


G.  SECTION  SUMMARY 


Vector  switch  decomposition  was  reviewed  from  the  viewpoint  of 
multi-rate  sampling.  For  particular  configurations,  it  was  demonstrated 
that  an  operational  algebra  could  be  used  in  lieu  of  switch  decompo¬ 
sition  and  thus  a  scalar  problem  can  retain  a  scalar  format.  It  was 
conjectured  that  all  open-loop  configurations  and  those  closed-loop 
systems  wherein  the  frame  times  are  ratioed  as  powers  of  2  constitute  a 
class  of  problems  on  which  the  operational  algebra  can  be  effectively 

applied.  The  algebra  requires  the  evaluation  of  nested  multi-rate 

r  T/4iT/2 

operations  such  as  [G3G1  J 


SECTION  IV 


A  USEFUL  ALGORITHM  FOR  THE  ANALYSIS  OF 
MULTI -RATE  CONTROLLERS 

A.  INTRODUCTION 

A  prevalent  trend  in  digital  control  systems  design  is  the  use  of 
several  computational  frametimes.  For  example,  both  the  F-18  and  Space 
Shuttle  have  three-rate  digital  controllers  (80,  40  and  20  Hz;  100,  50, 
25  Hz).  Moreover,  it  is  not  uncommon  to  find  sophisticated  simulations 
that  use  more  than  one  computer,  each  one  working  in  a  different  frame 
time  (e.g.,  see  the  case  study  of  Section  IX).  The  purpose  of  this 
section  is  to  present  an  algorithm  which  is  useful  for  the  analysis  of 
such  systems.  Starting  with  a  significant  generalization  of  the  "skip¬ 
sampling  theorem,"  we  first  comment  on  the  properties  of  the  solution 
and  from  these  deduce  an  algorithm  for  computing  two-rate  transfer  func¬ 
tions.  The  algorithm  requires  only  the  use  of  synthetic  division  as  the 
primary  analytical  tool.  Therefore  a  host  of  problems  associated  with 
alternative  methods  are  circumvented.  The  two-rate  algorithm,  when 
coupled  with  the  generalized  skip  sampling  theorem,  expedites  the  analy¬ 
sis  °f  multi-rate  systems. 

Examples  are  used  to  demonstrate  the  various  properties  of  the 
method  —  perhaps  the  most  important  being  that  the  algorithm  treats 
both  the  so-called  high-to-low  and  low-to-high  rate  cases,  as  limiting 
forms,  within  a  single  framework. 

B.  AN  IMPORTANT  IDENTITY 

Of  fundamental  concern  is  the  output  for  the  system  of  Figure  24. 
As  indicated,  it  is  permissible  to  insert  a  phantom  sampler  between  G 
and  the  T/M  sampler  if  T/NK  and  T/M  are  presumed  to  be  integer  related. 
A  slightly  less  general  version  of  Figure  24  wherein  K  -  M,  is  treated 
in  References  1  and  2.  Thus 

CT/M  „  [GrT/N]T/M  -  [GT/NKRT/N]T/M  (116) 
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G 


C 


G 


Figure  24.  A  Phantom  Sampler  Formulation  of  a 
T/N,  T/M  Sampling  Format 


where  M  and  N  are  rational  numbers.  In  Equation  116,  the  notation  of 
Reference  3  is  utilized.  Specifically,  G[’]  designates  the  z  transfer 
function  in  a  [•]  timeframe.  For  example. 


-  e-T/2 


Also,  knowledge  of  routine  operations  such  as 


(117) 


l  |T/2lT/6 


_  e-T/2 


is  presumed. 

For  the  special  case  M  ■  1,  K  -  1,  Figure  24  reduces  to  the  well- 
known  high-to-low  rate  transform  (for  example,  see  Reference  4). 

To  compute  the  entries  of  Equation  116  one  may: 

1)  Determine  R  by  table  lookup. 

2)  Replace  T  by  T/N. 

£ 

3)  Replace  z  (in  R)  with  z  . 

4)  Determine  G  by  table  lookup. 


*  <* .  ■ 


r  * 


5)  Replace  T  with  T/NK. 


6)  2  remains  the  same  (in  G). 

7)  Compute  CT^M  using  Equation  119. 


cI/“  ‘  nj  /r  G(p)  R<pK> 


2  dp 

z  _  pNK/M  p  p 


res 


G(p)R(pK) 

P 


(119) 
sT/NK 


In  Equation  119,  keep  in  mind  the  primary  requirements;  z  -  e 
NK/M  and  K  must  be  integers.  This  in  turn  requires  M  and  N  to  be 
rational  numbers.  An  example,  using  Figure  24,  will  clarify  the 
details . 

Example  1 

Let 


sT/6 


M  =  3 


2,  K 


s  +  2 


gt/6 


z  -  e“T/3 


s  +  1 


so  that 

cT/3 


rT/2  _ 


z3  -  e-T/2 


(z  -  e-T/3](z3  _  e-T/2) j 


,T/3 


p^z 


/p  (p  _  e-T/3)(p3  -  e-T/2) (2  -  p2)  P 

-L-  f  - _ 

211  J  Jy  (p-e“T/3) (p3_e~T/6) ^p2+  e-T/6p+e“2T/6)|)z_p2) 


dp 


res 


p-e-T/3  +  r6S 


+  res  |  .  +  res  ,,  T/,  - 

»e-T/6  p— e-T/6(l/2)+je-T/6(/372) 


p«e 


+  res 


p—e-T/6  (1  /2)-je-T/6  (/1 71) 


(120) 


The  task  is  to  evaluate  four  residues  in  the  z-plane,  even  though 
both  functions,  in  the  s-plane,  were  only  first  order.  Evaluating  the 
four  residues  and  placing  the  sum  of  the  four  terms  over  a  common  denom¬ 
inator  gives 


CT/3 


_ _ z^(z^  +  e~^/k)  _ 

(z  -  e~^T/3j(z  -  e“T/3j^z2  ■+■  e~T/3z  +  e~2T/3 j 


**  4  e"5T/6z3 


(z  -  e"2T/3)(23  -  e“T) 


Thus  one  starts  with 


cT/V/2 


{z  -  e-T/3)(z3  -  e-T/2) 


and  finishes  with 


cT/3  = 


z4  +  e-5T/6z3 


[z  -  e“2T/3)(z3  _  e-T) 


(121) 


z  =  esT/6  (122) 


z  =  esT/3  (123) 


leading  to  the  following  observations: 

1)  The  poles  of  Equation  123  can  be  directly  deter¬ 
mined  from  an  inspection  of  Equation  122  —  the 
root  locations  are  squared  since  the  ratio 
between  the  T/6  and  T/3  frame  times  is  2. 

2)  The  numerator  of  Equation  123  is  fourth  order  — 
one  may  adopt  the  viewpoint  that  the  entire 
analytical  effort  expended  on  the  evaluation  of 
Equation  119  has  the  net  effect  of  determining 
the  five  numerator  coefficients  of  Equation  123 
(NUM  =  z^  +  e~5T/6z3  +  qz2  +  Qz  +  Q). 


C.  OUTLINE  OF  A  NEW  PROCEDURE 


The  illustrative  example  suggests  an  alternative  to  evaluating  the 
residues : 

1)  Determine  the  denominator  polynomial  in  the  new 
time  frame  (T/M)  from  a  knowledge  of  the  pole 
locations  in  the  "old"  timeframe  (T/NK)  utilizing 
NK/M  as  the  ratio  between  the  two. 
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2)  Utilize  the  knowledge  of  the  poles  and  zeros  in 
the  T/NK  timeframe,  together  with  awareness  of 
the  pole  locations  in  the  T/M  timeframe,  to 
generate  the  zeros  in  the  T/M  timeframe. 

At  this  juncture  Step  1  is  clear  enough  but  a  definitive  method  for 
implementing  Step  2  is  needed.  We  shall  demonstrate  that  matching  the 
time  responses  of  Equations  122  and  123  furnishes  a  computationally 
efficient  method  for  implementing  Step  2.  Before  proceeding  to  this 
step,  the  reader  may  find  it  advantageous  to  survey  the  examples  of 
Appendix  A,  designed  to  provide  insight  into  the  setup  of  Equation  119 
using  Figure  24.  In  particular,  note  that  both  the  classical  "High-to- 
Low"  and  "Low-to-High"  rate  transforms  are  treated  as  limiting  cases. 

0.  IMPULSE  RESPONSE  MATCHING 


The  implementation  of  Step  2  above  by  matching  the  impulse  response 
is  described  next.  Given  an  nth  order  "transfer  function"  in  a  T/N  time 
frame,  the  transfer  function  in  a  T  time  frame  can  (computationally)  be 
found  by  matching  the  impulse  response.  Let 

NUM  -N  ,  -2N 
=  DEN  =  an  +  •••  +  an-iz  +  •••  +  an_2z 


+  •••  +  aQZ  +  • • •  ,  z  =  (124) 

where  n  is  the  order  of  the  system  and  N  defines  the  ratio  between  the 
time  frames. 

From  Equation  124,  form  a  vector  composed  of  n  +  1  terms  from  the 
successive  Nth  points  in  the  T/N  time  frame  transient  response: 

a  “  [an»  an-l »  *•*  ao] *  025) 

Postulate  a  transfer  function,  in  the  T  timeframe,  in  which  the  coeffi¬ 
cients  of  the  denominator  are  known  but  the  numerator  coefficients  are 


not.  Nothing  essential  is  lost  by  setting  the  denominator  lead  coeffi¬ 
cient  to  unity. 


GT  - 


-nzn  +  cn  -l2*1-1  +  •••  +  cp 
zn  +  +  •  •  •  +  bQ 


(126) 


A  vector  "b"  can  be  formed  from  the  known  coefficients  of  the  denomina- 


b  "  U»  bn-l»  bn-2»  *“»  b0 1* 


(127) 


Let  the  vector  "C"  represent  the  unknown  coefficients  of  the  numerator: 


[cn.  cn_i,  cn_2,  ***,  c0]' 


(128) 


The  solution  for  "c",  deduced  from  equating  the  partial  fraction 
expansion  of  Equation  126  to  the  correct  temporal  terms  of  Equation  124, 


an-l  an 


an-2  an-l 


32  ....  an 


(129) 


To  show  this,  observe  that  z-N  from  the  old  time  frame  is  z“*  in  the 
new  time  frame,  and  form 
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an  +  an-lz_1  +  an-2z_2  +  •••  +  a0z-n  + 


— 
sn  +  bn-izn_1  +  bn-2zn'2  +  •••  +  bQ  cnzn  +  cn_jzn-l  +  cn-2zn-2  +  •••  +  cQ 


an  anbn-l _ anbn-2 _ an  bQ 

cn-l“anbn-l  . 

an-l  . 


Equating  terms  gives 


cn  ~  an  =  ® 


cn  "  ar 

The  first  division  gives 

cn-l  "  anbn-l  =  an-l 


or 

cn— 1  *  an-l  +  anbn-l 

A  second  division  yields 


(130) 


(131) 


(132) 


cn-2  ~  anbn-2  ~  an-lbn-l 

cn-2 


an-2 

an-2  +  an-lbn-l  +  anbn-2 


(133) 


Continuing  for  n  divisions  yields  the  final  result,  which  can  be  placed 
in  the  matrix  form  of  Equation  129.  Several  examples  are  given  to  clar¬ 
ify  the  details. 


Example  2  (a  simple  check  case) 

Let 

C^/3  m  - -  1  +  ,9z~l  +  ,81z-2  +  .729z-3  +  •••  ,  z  ■  esT/3 

z  -  .9 

(134) 
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one  may  write 


c3 

c2 

C1 

_c0_ 


0 

4 

144 

1984 


0 

0 

4 

144 


U 

0 

0 

4 


0 

1 

0 

0 

-8 

4 

0 

-112 

112 

0 

-256 

384 

Summarizing  the  example,  it  is  seen 


CT/4 


or 


z2  —  4z  +  6 _ 

(z  -  2)(z2  -  2z  +  2) 


=>CT 
2»esT/ ^ 


4z2  +  112z  +  384 
(z  -  16)(z  +  4)2 


pT  4z  +  96 

(z  -  16)(z  +  4)  2,e-sT 


(139) 


We  conclude  the  examples  by  revisiting  Example  1. 


Example  5 


R  /  , 

1 

/  CT/V  cT/3 

1  T/2 

s+l 

s+2 

T/6  T/3 

Figure  25.  Example  5  (and  1)  Block  Diagram 


From  Figure  25, 


CT/3 


[/  1  \T/6/  !  \T/2 

T/3 

z  z^ 

~(z  -  e"2T/6)  (*3  .  e-T/2) 

T/3 
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or 


cT/6 


z4  _  e-T/3z3 

1  +  e-T/3z-l 


z^ 

-  e"T/2z  +  e"V6T 
+  e-2T/3z-2  +  ... 


cT/3 


C4Z^  +  C3z3  +  C2Z^  +  cjz  +  cq 
(z  -e"2T/3)(z3  _  e-T) 


(140) 


A  tedious  but  straight  forward  exercise  gives 


c4 

•P- 

o 

o 

o 

o 

_l 

1 

1 

c3 

a3  a4  0  0  0 

_e-2T/3 

0 

c2 

“ 

a2  33  a4  0  0 

0 

— 

e-5T/ 6 

ci 

al  a2  a3  a4  0 

i 

n> 

I 

H 

0 

_c°- 

a0  al  a2  a3  a4 

e-5T/3 

0 

where 

a  »  [a4,  a3,  a2,  aj  aj]' 

«  [1,  e-2T/3f  e-4T/3  +  e~5T/6f  e-T  +  e-2T  +  e-31/2. 


e-5T/3  +  e-8T/3  +  e-13T/6]' 


(141) 


(142) 


E.  SECTION  SUMMARY 

The  system  of  Figure  24  can  be  evaluated  using  Equation  119,  which 
Is,  in  itself,  a  significant  generalization  of  the  "skip-sampling 
theorem."  In  particular,  M  and  N  need  only  be  rational  numbers.  This 
relaxed  requirement  on  M  and  N  permits  the  treatment,  as  limiting  cases, 
of  the  so  called  High-to-Low  and  Low-to-High  rate  transforms  within  the 
same  framework.  These  points  are  called  out  in  Appendix  A. 
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Once  set  up  via  Equation  119,  the  problem  reduces  to  the  evaluation 


of 


CT/M  „  [cT/NK]T/M 

One  may  avoid  using  residue  theory  by  doing  a  powers  series  expansion  of 

^T/NK  an(j  equating  the  NK/M  spaced  points  in  the  T/NK  timeframe  to  the 

appropriate  temporal  points  in  the  T/M  time  frame*  This  leads  to  a 

T/M 

closed-form  solution  for  the  C  numerator  coefficients.  In  this 

T  /M 

regard,  .  it  is  important  to  appreciate  that  the  poles  of  C  '  can  be 

T  /Wlf 

obtained  from  the  known  poles  of  CL/  "by  inspection." 

The  combination  of  the  generalized  skip-sampling  algorithm  and 
impulse  response  matching  eases  significantly  the  computational  burden 
encountered  in  the  analysis  of  multi- rate  systems. 
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SECTION  V 


HOLT I -NATE  FREQUENCY  RESPONSE:  SWITCH  DECOMPOSITION 
CONTRASTED  WITH  SCALAR  ANALYSIS 


A.  INTRODUCTION 

Vector  switch  decomposition,  discussed  in  Section  III,  provides  a 
general  framework  for  the  frequency  response  of  multi-rate  systems  (Ref¬ 
erences  1,  2).  However,  the  dimensionality  problems  associated  with 

switch  decomposition  encourage  the  development  of  a  separate  theory  for 
the  class  of  systems  discussed  in  Section  III.  We  first  review  the  gen¬ 
erally  applicable  multi-rate  frequency  response  (as  presented  in  Refer¬ 
ence  2)  in  order  to  make  the  additional  dimensionality  problems  clear. 
Following  this,  a  less  general  method  is  developed.  This  "scalar" 
method  is  applicable  to  the  class  of  problems  discussed  in  Subsection  E 
of  Section  III  and  will  be  applied  to  the  simulation  case  studies  of 
Sections  VIII  and  IX. 

B.  REVIEW  OP  MULT I -RATE  FREQUENCY  RESPONSE, 

SWITCH  DECOMPOSITION  MODEL 

Let  the  general  multi-rate/multiple-order  open-loop  system  of  Fig¬ 
ure  26  have  a  sine  wave  input.  In  Figure  26, 


C1  - 

(GiR<*)e 

(143) 

CT/N  - 

[G(GlRa}6] 

(144) 

where  a,  6  represent  sampling  schemes  with  a  basic  period  of  T  seconds. 
For  example,  a  might  represent  a  multiple-order  sampling  format;  0  might 
represent  a  multi-rate  and/or  pseudo  measurement  format  (see  Refer¬ 
ence  2 ) . 


Figure  26.  Multi-Rate/Multiple-Order 
Open-Loop  System 


Using  switch  decomposition,  Figure  26  takes  on  the  representation  of 
Figure  27.  Clearly, 

cj  -  W2(w2*G1W1)T(wUr)T  (145) 

and 


CT/N  -  (gW2)T/N  (w2*GiW1)T  (wur)T  (146) 


Figure  27.  Open-Loop  System  with  Switch  Decomposition 


If  a  represents  multiple-order  sampling  and  6  a  pseudo  measurement 
format  using  multi-rate  sampling,  the  switch  decomposition  modeling  com¬ 
ponents  might  appear  as : 
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wx 


[1. 


a-T0S 


-Tis 


-T2s 


(147) 


and 


wl* 


eTos 

eTls 

eT2s 


(148) 


W2  = 


W2 


1  0  0 

o  e-sT/3  0 

0  o  e-2sT/3 


(149) 


w2* 


1 

esT/3 

e2sT/3 


(150) 


Comparing  Equacion  146  with  Equation  13,  we  see  that  only  one  new 
facet  has  entered  the  problem,  namely  (Wj  R)^  replaces  R^.  Consider, 
therefore,  a  generic  component  of  (Wj.^r]T  —  for  instance  the  scalar 
(0ATsr)T,  where  0  <  A  <  1.0. 

For  R  *  sin  bt  and  A  zero  (e^s  =  1),  the  output  equation  (Equa¬ 
tion  146)  becomes 


CT/N 


(GW2)T/N(W2*G1W1  )T 


_ zN  sin  bT _ 

Z2N  _  (2  C08  bT)zN  +  1 


z  -  esT/N 

(151) 


T 

where  R  is  described  in  terms  of  a  N/T  samples  per  second  model.  For 
the  sake  of  brevity  write  Equation  151  as 
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Expand  the  right-hand  siae  of  Equation  152  in  partial  fractions: 


CT/N 


N-l 


£ 

n«0 


sin  oin(T/N)  +  Bnz[z  -  cos  mn(T/N)] 
z^  -  [2  cos  cun(T/N)]z  +  1 


For  non-zero  A,  we 
Equation  153  becomes 


+  [ 


Terms  due  to  modes  of  G 


T/N„T] 
A  SJ 


use  the  advanced  z-transform  on 


<.1TsR) 


(153) 

r 

and 


CT/N 


T/N  T  zN{(sin  bAT)zN  +  sin  [b(l  -  A)T] 
A  B  z2N  .  (2  cos  bT)ZN  +  1 


(154) 


T  /N  T 

Assume  that  responses  in  the  modes  of  G,  G„  approach  zero  as  t  +  °°, 

A  D 

i.e.,  that  all  modes  are  stable.  In  Equations  153  and  154 

2nn 

oijj  *  b  +  ,  n  “  0, 1,  2 . N-l 

(or  one  may  prefer  the  n  “  0,  etc.,  definitions  of  Section  II).  The 

steady-state  waveform,  at  the  sampling  instants,  can  be  written  as 


C(t) 


T/N 


N-l 

£ 

n-0 


T/N 


(An  sin  u>nt  +  Bn  cos  u>nt) 


(155) 


At  this  point,  pursue  the  analysis  of  Section  II  keeping  in  mind  that 
the  only  difference  resides  in  the  use  of  an  advanced  z-transform  for 
the  sine  wave  input.  The  analysis  proceeds  down  to  the  equivalent  point 
of  Equation  152  where  a  division  by  z  sin  wk(T/N)  and  replacement  of  k 


with  n  occurs 
limiting  form: 


At  this  point  we  pick  up  the  development  and  find  the 


.  ~T/N_T 

A„+JBn  -  CA  Gg 


zN~1[(sln  bAT)zN  +  sin  b(l-A)T)  *  ~t2  cos  mn(T/N)]z+l 


z-U^(T/N) 


sin  «ip(T/N) 


z2N-(2  cos  bT)zN+l 


Evaluating  Equation  156  at  z  -  14-Gin(T/N)  produces 


An+JBn  "  CA  CB  2“  lgii»n  (T/N ) 


{sin  bAT(cos  WpT  +  j  sin  wnT)  +  sin  b(l  -  A)T} 

_ «  [cos  Up(T/N)  +  j  sin  uip(T/N)  -  cos  oy,  (T/N ) ) 

[sin  uin  (T/N )  JN  (cos  nipT  +  j  sin  u>nT  -  cos  bT] 


Simplifying  Equation  157  gives 


An+JBn 


rT/NrT 
gA  gB 


z-l*w„(T/N) 


(sin  bAT)(cos  w^T)  +  sin  b(l  -  *  >T  +  J(sin  wnT)(3ln  bAT) 
N  a  In  b 


-  G 


T/N-T 
A  °B 


z-ltUp  (T/N) 


sin  bT  cos  bAT  +  J  sin  wnT  sin  bAT 
N  sin  bT 


T/N.T  ,  cos  bAT  -t-  1  sin  bAT 

A  gB  2-i4uin(T/N)  N 


-  G 


T/NrT I  .  eJbAT 

'A  ’Bfz-l»un(T/N)  N 


T/N  T 

Ga  gb 


*-l*«*n  (T/N) 


e^U-jb 


-Us.n(T/N) 

Cl  56) 


(157) 
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Thus,  the  only  new  element  added  is  the  factor  eAS  evaluated  at  the 

input  frequency-  Since  A  is  generic,  we  draw  the  same  conclusion  for 

T 

every  other  element  of  (W^R)  and  hence  arrive  at  the  final  steady 
state  result: 

,  T 

n-1 

CT/N  „  (GW2^T^N(w2*Glwl)T(wl*R)T  “  S  (An  sin  “nt  +  Bn  cos  wnt) 

n*0 


where 


An  +  JBn 


(GW2)T/N|zMT/N 


i  z=es^ 


(w2*GlwlJT  WU|  .. 

z  =  l5.un(T/N)  *  gp^^n1  *s=jb 

z  =  l*  bT 


(159) 


“n  =  b  + 


n  “  0,1 . N-1 


(160) 


Letting  N  ■*  °°  gives  the  "continuous"  result 


GW2|  _  z=e& A 

An  +  JBn  *  —  (w2*Glwl)T  w1*  ..  n=0»  >  *2>  *•* 

ls=jwn  |r1^u>nT  s“jb 

z=l4bT 


(161) 


T 

In  Equations  159  and  161,  an  option  exists  with  regard  to  C W2*G1W1 ^  as 

it  may  be  evaluated  at  either  l$.bT  or  l$.u)nT.  However,  there  is  no 

option  with  regard  to  since  it  must  be  evaluated  at  s  =  jb.  Clear- 

^  T 

ly,  the  most  efficient  procedure  is  to  evaluate  (w2*^lwl'  at  z  =  l&bT 
and  at  s  ■  jb,  save  the  result  in  a  polar  format,  and  then  evaluate 
(GW2/N)^N  —  or  GW2/T  —  over  the  appropriate  range  of  n. 

It  can  be  appreciated  that  the  results  of  Section  II  are  contained 
in  Equations  159  and  161,  since  Section  II  considers  only  the  special 
case  of  Wj  *  1.0  (a  scalar).  Also,  note  that  Wj  and  W2  are  symbolic 
representations  that  can  represent  either  multi-rate  or  multiple-order 
sampling  formats  (see  for  example.  Reference  2).  However,  in  arriving 


at  the  spectrum  of  the  continuous  variable,  uniformly  spaced  time 
samples  in  a  T/N  time  frame  were  used. 


C.  A  CLOSED-LOOP  APPLICATION 

The  ability  to  cast  any  multi-rate  or  multiple-order  sampling 

sequence  into  a  switch  decomposition  format  makes  the  analysis  of 

closed-loop  systems  a  straightforward  task.  Consider  the  system  of  Fig- 

T 

ure  28.  Proceed  as  in  the  single-rate  case.  Since  C  ■  G^W^E  »  one  must 
T 

first  solve  for  E  : 

ET  -  (wUr)T  -  (wUG2W2)T(w2aG1W1)T  ET  (162) 


or 


eT  -  [i  +  (wuG2W2)T(w2jtG1W1JT]-1(wuR]T  (163) 


Figure  28.  Closed-Loop  System  with  Switch  Decomposition 


Therefore 

CT/N  -  (giWi)T/N[i  +  (wuG2W2)T(w2jkG1W1)T]-1  (vur)T 


(164) 


73 


The  coefficients  for  the  steady-state  waveform  are  then 


^n  +  JBn 


G1W1 


s-jo)n 


[I  +  (wuG2W2)T(w2ikG1W1)T] 


-1 


z=e 


sT 


W, 


z-l4.bT 
or  z=14u)nT 

°»  +1*  ±2, 


s-jb 

(165) 


Section  IV  of  Reference  2  gives  a  detailed  example  on  the  interpretation 

of  Equation  165  and  we  will  not  repeat  it  here.  However,  an  important 

point  made  in  Reference  2  will  be  re-emphasized  here  —  namely  the 

presence  of  Wj^l  forces  the  use  of  "multiple  Bode  plots",  one  for 

I  s=  ib 

each  component  of  »!*•  This  is  yet  another  facet  of  the  dimensionality 
problems  associated  with  switch  decomposition. 

For  a  particular  class  of  problem  (discussed  in  Section  III),  the 
dimensionality  problems  can  be  avoided  by  using  a  more  direct  approach 
which  follows. 

D.  A  DIRECT  APPROACH 

As  pointed  out  in  Section  III,  open-loop  systems  and  a  class  of 
closed-loop  multi-rate  systems  can  be  analyzed  without  recourse  to 
switch  decomposition.  Such  is  the  case  with  the  generic  system  of  Fig¬ 
ure  29. 


rT/M 

ft  , 

/  * 

n  ~ 

T/N 

T/K 

t>  2 

T/M 

Figure  29.  A  Three-Rate  Scalar  Configuration 


Clearly, 


CT/M 


{g2[g1rT/n]t/k} 


T/M 


(166) 
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It  would  be  fortuitous  Indeed  for  the  ratios  of  N,  K,  M  to  be  such 
that  all  the  operators  in  Equation  166  "operated"  through.  However, 
knowing  that  N,  K,  and  M  are  rational,  one  may  use  the  algorithm  of  the 
previous  section  to  evaluate 

gA/K  “  [G1RT/N]T/K  (167) 


The  poles  of  RT^N,  In  a  T/N  time  frame,  will  survive  Intact  In  the  T/K 
time  frame,  indicating  that  the  input  R  can  be  "tracked"  through  the  T/N 
and  T/K  operations.  In  a  like  manner 


CT/M 


(02c I/K) 


T/M 


(168) 


can  be  evaluated  and  the  poles  of  R,  in  a  T/M  time  frame,  can  be  clearly 
identified.  This  suggests  a  possible  method  for  avoiding  the  switch 
decomposition  frequency  response  by  factoring  out  a  sine  wave  in  the  T/M 
time  frame,  thus  placing  the  investigator  in  a  position  to  apply  single¬ 
rate  frequency  response  results  directly.  This  viewpoint  is  flawed,  as 
can  be  appreciated  by  studying  a  particular  case.  Suppose  K/N  is  an 
integer  and  M/K  is  an  integer.  Equation  166  becomes 

CT/M  »  c^GfV/N  (169) 

and  clearly  the  poles  of  R  appear,  in  the  T/M  time  frame  (z  ■  esT^M)  in 
terms  of  zM^N,  i.e., 

z2  -  2  cos  (bT/N)z  +  1  z2M/N  -  2  cos  (bT/N  )zM^  +  1  (170) 

Thus,  the  frequency  response  of  requires  M/N  sinusoids  in  order  to 
match  the  output  samples  in  a  T/M  time  frame.  That  is,  one  sine  wave, 
in  a  T/M  time  frame,  is  not  a  sufficient  descriptor  of  the  sampled  T/M 
steady  state  waveform. 

One  may  obtain  a  consistent  evaluation  for  the  spectral  content  of 
the  output,  in  the  context  of  generating  the  minimum  set  of  sinusoids 
which  will  exactly  match  the  output  samples,  by  careful  use  of  the 
generalized  skip-sampling  theorem  developed  in  Section  IV.  We  will  use 
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examples  to  demonstrate  the  methodology  and  proceed  to  discuss  a  defini¬ 
tive  set. 

Example  #1; 


Figure  30.  Illustrative  Example  #1 

T/nT/Ni1 


CT  .  [GlRT/N]r  -  [g^7  R7  ] 


(171) 


For  ease  of  identification,  let  Gj  =  (1/s  +1),  R  *  sin  bt. 
Therefore, 


CT 


z  sin (bT/N) 


-  e“T/N  z2  -  2  cos (bT/N )z  +  1 


(172) 


— T /N 

Since  the  ratio  between  the  frame  times  is  N,  the  root  at  z  =  e  in  a 

— T 

T/N  time  frame  goes  to  z  *  e  in  a  T  time  frame.  The  roots  correspond¬ 
ing  to  the  sine  wave  in  T/N,  where 


z  =  cos  (bT/N)  +  j  sin(bT/N)  -  l$.bT/N 


jump  to  zN  »  (lzJ_bT/N)N  or  z  -  14-bT.  Thus  CT  will  have  the  form 


(173) 


CT  - 


N(z) 


(z  -  e_,r)[z2  -  (2  cos  bT)z  +  l] 


z-esT  (174) 


and  it  can  be  appreciated  that  only  one  sine  wave  is  needed  to  match  the 
steady  state  response  (only  because  of  fast  input-slow  output  sampling 
format) • 


Figure  31 •  Fast  Output,  Slow  Input  Example 


CT/N  _  [G  RT]T/N 


gT/NrT 


_ z _ z^sin  bT _ 

(z  -  e-T/N]  [Z2N  -  (2  COS  bT) z^  +  l] 


(175) 

z.esT/N 


It  is  apparent  that  the  sum  of  N  sinusoids  is  required  to  exactly  match 
the  T/N  output  samples. 


Example  #3: 


Figure  32.  Three  Rate  Example 


Again,  for  ease  of  identification,  let  -  l/(s  +  1)  and  G2  «  1/ 

(s+2). 

T/3 


T/2, 


cT/3  -  {g2[g1rt/16]  } 


(176) 


The  inner  operation  will  yield,  using  the  method  of  Section  IV, 

t /2 

z _ z^  sin(bT/16) _ 

z  -  e-T/16  z2  _  j2  co8(bT/16)]z  +  1 


U 


_ MU _ 

e-T/2) [z2  _  [2  cos (bT/2)] z  +  1 


_ sT/2 


(177) 


since  the  ratio  between  the  frame  times  is  eight.  The  next  operation 
must  be  initiated  in  a  T/6  time  frame,  which  gives  a  form  like 


Gi/6[GlRT/16]T/2 


N  (z  3 ) 


(z3  -  e“T/2)(z6  -  [2  cos (bT/2 ) ] z  +  1 


T/3 

(178) 


;  ratio 

between 

T/6 

and 

T/3 

is  two 

,  and 

looking  at  the  T/6  roots: 

2  *» 

-2T/6 

e 

» 

2 

z 

- 

-2T/3 

e 

(z  -  e_2T/3)  in  T/3 

3 

2  « 

-T/2 

e 

» 

6 

z 

- 

-T 

e 

(z3  -  e“T)  in  T/3 

(179) 

3 

z  = 

14-  bT/2 

> 

6 

z 

1  bT 

[z6  -  (2  cos  bT)z3  +  1]  in 

T/3 

CT/3  has  a  form 

cT/3  „  - N.l(z) _ 

(z  -  e-2T/3)(z3  ~  e_T)(z^  -  (2  cos  bT)z3  +  l) 

Thus  the  "fragmentation"  of  the  sine  wave  through  T/16,  T/2,  and  T/3 
time  frames  requires  the  use  of  three  sinusoids  to  match  the  output 
steady  state  waveform  (at  the  sampling  instant). 

In  this  example,  one  can  therefore  use  single  rate  results  (avoid 
switch  decomposition)  by  setting  up  a  new  transfer  function  via 


(180) 

-esT/3 
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_ Nt  (z) _ 

(z3  sin  bT](z  -  e~2T/3j(23  -  e-^) 


z6  -  (2  cos  bT)z3  -  1 


(181) 


Gi(z) 


sin  bT 


-  (2 


(182) 


z  sin  bT 


sesT/3  z2  _  (2 

cos  bT)z  +  1 


(183) 


I  z=esT 


An  +  j»n  "  i  Gl<z> 


UT/3 


,  ton  =  b  +  (2irn/T) 

n  =  1,  +1,  or  0,  1, 2 


Clearly,  the  same  problem  which  surfaced  in  the  switch  decomposition 

frequency  response  also  surfaces  here.  Namely,  in  switch  decomposition 

w*  was  evaluated  at  the  input  frequency  b;  it  was  not  permissible  to  let 

T/3 

its  value  range  with  wn  In  Gj  ,  it  is  necessary  to  keep  the  scale 
factor,  sin  bT  in  the  denominator  of  G^,  fixed  at  the  value  b. 

These  examples  make  it  clear  that  the  use  of  this  scalar  theory,  as 
an  alternative  to  switch  decomposition,  requires  insight  into  the  manner 
in  which  th ;  poles  are  manipulated  between  successive  time  frames. 

E.  SECTION  SUMMARY 


The  frequency  response  of  multi-rate  systems,  using  a  switch  decom¬ 
position  format  was  reviewed  with  a  sinusoidal  input.  Decomposition 
of  the  input  (w*r)^  requires  a  number  of  "Bode  plots"  equal  to  the 
dimension  of  W^.  An  alternative  "scalar  approach,"  using  the  tool 
developed  in  Section  IV,  was  presented.  This  approach,  where  applic¬ 
able,  is  dimensionally  more  attractive  than  vector  switch  decomposition. 
It  does  require  a  keen  awareness  (on  the  user's  part)  of  the  number  of 
sinusoids  required  to  match  the  sampled  output  steady-state  waveform. 


SECTION  VT 


SIMULATION  MODELS 


A.  INTRODUCTION 

There  is  a  natural  extension  of  the  previous  results  which  make 
switch  decomposition  (or  its  equivalent)  and  frequency  response  tech¬ 
niques  ideal  for  the  error  analysis  of  simulations.  Before  discussing 
that  topic  in  explicit  detail,  we  will  first  develop  models  for  treating 
a  given  computer  code  under  the  assumption  that  it  introduces  an  inher¬ 
ent  throughput  delay.  The  limited  ability  of  "time  advanced"  digital 
filters  to  compensate  for  throughput  delay  is  briefly  treated.  A  key 
idea,  the  use  of  a  zero-order  hold  to  model  a  buffer  storage  register 
between  two  computers  (working  in  different  frame  times)  and/or  differ¬ 
ent  elements  of  a  computer  program,  is  introduced  first. 

B.  MODELING  A  BUFFER  REGISTER 

The  use  of  a  ZOH  to  model  a  storage  register  is  depicted  in  Fig¬ 
ure  33.  There  the  output  of  the  filter  G(z)  is  stored  in  an  intermedi¬ 
ate  buffer  register,  Mq.  The  output  C  is  given  as 

C  -  M0GTRT  (184) 


and  clearly 


Figure  33«  Model  of  a  Buffer  Register 


The  utility  of  this  model,  therefore,  lies  not  in  the  situation  where  C 
is  sampled  at  (1/T)  Hz,  but  rather  when  some  other  type  of  output  opera¬ 
tion  is  performed.  For  example,  one  may  visualize  C  as  being  the  input 
for  the  next  (serial)  recursion  operation  to  be  performed  in  a  different 
frame  time,  say  T/3  (refer  to  Figure  34).  Express  C  as 


C  -  M3G2/3  [MjG^R1] 


T/3 


M3G!/3Ml/3Gy 


(187) 


where  the  subscripts  on  the  data  holds  are  used  to  indicate  the  integer 
values  of  their  frame  time,  i.e., 


M1 


s 


m3 


(138) 


To  illustrate  the  use  of  Figure  32,  suppose  a  recursion  equation  like 


Xk  "  -*2X^.1  +  • 4Rfc  +  *3Rk-l 


(189) 


in  a  T  frame  time  is  followed  by  another  recursion  equation  in  a  T/3 
time  frame: 


Figure  34-  A  Particular  Serial  Operation 

The  transfer  functions,  with  z  commensurate  with  the  indicated  time 
frames,  are: 


_T  . 4z  +  .3  -T/3 

G1  =  z  +  .2  *  °2  = 

Equation  187  in  terms  of  a  T/3  time  frame  is: 

c  ,  1^*22  . 

s 


.  5z 

z  -  .  6 


(191) 


R(zJ) 


— -gl— — .  —5 _ JL^sCii  /». 4z.lt  ,i\R, 

s  z  -  .6  s  \  z3  +  .2  |R(Z 

_ jT.f-1/.3  .  — 5 - |z2_+  z_±  ,l|  /..4z3.+  .,3\  3 

s  2  ~  '6|  z2  I  z3  +  .2  / 


(192) 


(193) 


2  2 

The  term  (z  +  z  +  IJ/z  models  the  fact  that  recursion  Equation  19 0 
asks  recursion  Equation  189  for  an  input  every  T/3  seconds,  but  picks  up 
the  same  number  two  out  of  every  three  times. 

C.  MODELING  THROUGHPUT  DELAY 

Multiplication  and  addition  of  a  recursion  equation  require  a  finite 
amount  of  machine  time  for  their  execution.  Thus,  it  is  necessary  to 
understand  the  timing  of  a  simulation.  For  example,  the  recursion  equa¬ 
tion  of  Equation  194  can  be  machine  executed  in  several  different  ways: 


Yk  "  a  l^k-l  +  a2Yk-2  +  a3?k-3  +  boxk  +  blxk-l 

a)  Bring  in  the  new  x,  ,  do  five  multiplies  and  four 

adds,  and  (as  quickly  as  possible)  output  the  new 

value  of  y^.  There  will  be  a  throughput  delay, 

say  T  seconds, 
o 


(194) 


TK 


b)  Same  as  a)  except  y^  is  not  updated  until  the 
start  of  the  next  machine  cycle.  There  will  be  a 
throughput  delay  equal  to  the  frame  time,  T  sec¬ 
onds* 

Yet  another  alternative  is  to  perform  all  the  multiplies  and  adds  asso¬ 
ciated  with  the  past  values  as  a  "background"  computation.  Then,  at  the 
start  of  the  new  frame  time,  one  need  only  bring  in  the  new  value  of  xk> 
perform  one  multiply,  add  using  the  "background"  number,  and  output  y^ 
as  quickly  as  possible.  This  would  minimize  throughput  delay.  Explic¬ 
itly, 


Yk  =  boxk  +  bk 


(195) 


bk 


alyk-l  +  a2yk-2  +  a3yk-3  +  blxk-l 


(196) 


Thus,  while  the  delay  in  updating  x^  and  outputting  y^  as  quickly  as 
possible  must  be  minimized,  say  in  T^  seconds,  the  background  can  be 
done  in  a  more  leisurely  manner  provided  it  can  be  "fitted  in"  in  the 
remaining  T  -  T^  seconds. 

A  model  which  fits  all  three  cases  is  shown  in  Figure  35.  There  the 
nonsynchronous  sampler  model  of  Section  III  is  used  to  model  the 
throughput  delay.  This  model  envisages  that  the  delay  inherent  in 
either  the  time  required  to  multiply  and  add  (or  perhaps  a  deliberate 
delay  of  one  frame  time)  can  be  lumped  as  a  skewed  data  operation  after 
the  data  hold.  Thus  the  recursion  modeled  by  Gj  (z)  is  done  instantan¬ 
eously  and  passed  to  a  storage  register,  and  its  output  is  out  of  step 
with  the  master  clock  by  T0  seconds.  This  is  depicted  in  Figure  36. 

A  final  comment  is  in  order.  When  an  element  of  a  simulation  is 
modeled  as  a  recursion  equation  (computer  code),  a  single,  distinct 
frame  time  is  implied.  This  means  that  the  pulse  transfer  function  used 
to  model  it  must  have  the  same  frame  time  operation  on  both  input  and 
output. 
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Figure  35.  Buffer  Register  Model  with  Throughput  Delay 


xT 

T  + 

nT 

T  +  nT  +  T0 

eT 

Figure  36.  Timing  Diagram 

D.  SWITCH  DECOMPOSITION  MODEL 

In  Equation  187,  it  was  possible  to  operate  through  the  T  time  frame 
with  the  T/3  operator.  This  will  not  always  be  possible  and  we  may, 
once  again,  have  to  resort  to  switch  decomposition. 

With  switch  decomposition,  one  is  in  a  position  to  evaluate  configu¬ 
rations  such  as  the  one  shown  in  Figure  37  (the  vector  switch  decomposi¬ 
tion  is  shown  in  Figure  38)  • 

From  Figure  38, 

C  -  M2W2(W2*G2W2)T  (w2*M3W3}T(W3*G1W3)T  (w3*r)T  (197) 

The  evaluation  of  terms  such  as  (w2aM3W3)^  is  routine;  our  concern  is 
with  terms  like  (w3^G^W3)T  when  gT^  is  given  as  a  software  specifics- 
tion  (e.g.,  computer  algorithm).  Focusing  on  (W^  G^W^)  »  we  find 
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Figure  37.  Example  Two-Rate  Open-Loop  System 


Figure  38.  Switch  Decomposition  Formulation  for  Figure  35 


(w3*GiW3)t 


—  - 

\ 

1 

I 

esT/3 

gJi  |  e“8T/3  j  e-2sT/3n> 

e2sT/3 

1 
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1 

e-sT/3 
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01 
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e2sT/3 

eaT/3 

■  I 

(198) 
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T 


(W3*GlW3): 


Gl 


l  /  e-sT ,^2sT/3;  e-sT'eaT/3S> 

"'V 

(&IV  i  \  .-t/.2.t/3; 


/J2sT/3)  fefT/3' 


Jk—  T 
I 

1 ; 


(199) 


Observe  the  super-diagonal  terms*  They  are  simply  appropriate  sub- 

-sT 

diagonal  terms  multiplied  by  e  .It  is  only  necessary  to  compute  the 

T 

first  column  entries  in  order  to  completely  determine  (W3  ciw3)  • 

T/3 

Next,  observe  that  the  problem  statement  stipulates  ,  perhaps 

3  T/  3 

defined  through  a  substitution  rule*  If  e  is  defined  as  z  we  have: 


(w3*GiW3]t 


I 

1 


rT/3 

G1 

z-3(z2(.I/3) 

z-3(zG’f/3j 

zGT/3 

cT/3 

z-3(z2GT/3) 

2rT/3 

G1 

rT/3 

G1 

2  a  gST/ 3 

(200) 


Thus  one  may,  for  all  the  super  diagonal  terras,  factor  out  1/z  (in  a  T 
time  frame).  Although  each  terra  in  Equation  200  can  be  evaluated,  using 
the  algorithm  of  Section  IV,  only  the  first  column  entries  need  be  com¬ 
puted  (i.e«,  the  three  entries  of  column  one  are  sufficient  to  define 
all  nine  entries)  . 

E.  TINE-AOVAlfCED  DIGITAL  FILTERS 

In  synthesizing  a  recursion  equation  for  a  particular  element  of  a 
simulation,  the  designer  can  often  compensate  for  part  of  the  throughput 
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delay  introduced  by  the  digital  computer.  To  see  this,  consider  the 
task  of  implementing  a  low-pass  filter  via  a  recursion  equation  (refer 
to  Figure  39).  One  way  to  do  it  ia  to  visualize  the  input  to  the 
(sampled)  filter  as  being  smoothed  with  a  data  hold  (of  the  designer's 
choosing)  and  then  advanced  in  time  by  the  operator  e^8^.  If 
G  •  a/(s  +  a)  and  M  is  a  ZOH,  compute  the  pulse  transfer  function  as 


T 

G(z)  - 


I  -  s-T  .a.!/ — a— \ 

s  Is  +  a/ 


(201) 


or 


jeAsT 

1  _  _ 1 _ 

1 

s  s  +  a 

(l  -  e~AaT)z  +  fe-AaT  -  e~aT) 


_  e-aT 


(202) 


z  -  e 


It  can  be  appreciated  that  a  full  time  period  advance,  which  is  not 
physically  realizable  in  the  analog  world,  can  be  realized  on  the  dis¬ 
crete  domain,  since 


Cn  *  e~aTCn-i  +  ( 1  “  e_aT)^n 


(203) 


Figure  39*  A  "Time  Advanced"  Filter  Section 


f 


can  be  realistically  approximated  by  a  very  fast  computer.  This  method 
works  well  for  low-pass  sections  but  there  is  little  to  be  gained  by 
trying  it  on  washout  networks  or  high  passes.  To  see  this,  consider  a 
limiting  form,  such  as  a  differentiator: 


G(z) 


s 


eAsT 


(204) 


Here  any  A  *  0  results  in  a  physically  unrealizable  component.  Of 
course,  one  may  utilize  a  higher-order  data  hold,  for  instance,  the 
"slewer ." 


G  (z) 


(205) 


When  A  »  0,  the  smoother  is  called  a  slewer;  when  A  ■  1  the  smoother  is 
recognized  as  the  triangular  data  hold.  Clearly, 


G(z) 


1  fz  -  112  eAsT  1 
T  z2  s 


(206) 


regardless  of  the  value  of  A. 

We  close  this  facet  of  the  discussion  with  a  look  at  an  Integrator 


when  the  smoother  is  a  time-advanced  slewer  or  a  time-advanced  ZOH  and 
compare  the  results  with  the  Tustin  transforms  of  1/s. 


S lever : 


G(2)  -  ^ 

Tz2  s3 


(2  -  1)2  #  2  [(a2/2)  +  ( 1/2  -  A2  +  a)z  +  ( 1/2  -  A  +  A2) 
Tz2  (z  -  l)3 


(207) 


For  a  full  period  advance 


G(z) 


I  (Hi) 


(208) 


ZOH: 


G(z)  - 


z  -  1 

[esT 

T 

2-1 

z 

s2] 

m 

z 

(«  -  0: 


(209) 


For  a  half  period  advance  (unrealizable  if  A  ■  1) 


G(z) 


TuBtln  Transform: 


G(z) 


Tus tin 


1  (z  +  1) 

2  z  -  1 


(210) 


(211) 


Thus  Tustln  transform  of  1/s  Is  the  same  as  smoothing  with  a  ZOH, 
advanced  in  time  by  half  a  period,  or  a  slever  with  an  advance  of  T. 

F.  DIFFERENCE  EQUATIONS  FOR  MODELING 
AN  INTEGRATOR 

A  frequently  used  technique  for  deriving  recursion  equations  from  a 
given  transfer  function  is  to  rearrange  the  transfer  function  in  terms 
of  "nested"  values  of  1/s .  For  example,  the  network 


can  be  written  as 


G(s) 


a2/s  +  aj/s2  +  a0/s3 
1  +  b2/s  +  b|/s2  +  b0/s^ 


[a0/s  +  aj ] 1/s  +  a2 
(b0/s  +  b^)l/s  +  b2  j 1 / s  +  1 


(213) 


This  approach  emphasizes  the  importance  of  having  a  viable  discrete 
representation  of  the  integrator.  Indeed,  this  has  been  the  prime  focus 
of  the  classical  integrator  approach.  That  is,  given  an  input  x,  the 
integrator  produces  x  and  one  may  tabulate  a  set  of  classical  integra¬ 
tion  algorithms,  all  of  which  can  be  derived  using  the  approach  of  the 
previous  section  (see  Table  2).  In  the  remainder  of  this  report,  it 
will  be  assumed  that  the  recursion  equations  are  given  and  the  focus  is 
on  evaluating  of  the  fidelity  of  the  overall  simulation. 

G.  SECTION  SUMMARY 

Analytical  models  for  incorporating  recursion  equations  into  a 
framework  which  accounts  for  throughput  delay  were  developed.  The 
switch  decomposition  model  of  a  simulation  element  was  shown  to  have  a 
recursive  pattern  which  required  the  evaluation  of  only  the  first  column 
of  the  describing  matrix.  From  this  column  vector  all  other  entries  of 
the  matrix  can  be  generated.  A  fundamental  step  was  the  use  of  a  ZOH  to 
model  a  buffer  storage  register.  The  use  of  time  advanced  digital  fil¬ 
ters  to  compensate  for  throughput  delay  was  briefly  discussed  and  the 
classical  equations  for  integrators  were  summarized  in  tabular  form. 


TABLE  2.  CLASSICAL  INTEGRATORS 


INTEGRATION  EQUATION 

NAME 

*r  “  Xn_1  +  T(2Xn_1  -  Xn) 

— 

Sr  *  (3^1  -  Xr) 

— 

Sr  "  ^  +  T^j) 

Euler 

Sr  “  Sn_1  +|  (X„  + 

Trapezoidal 

Sr  -  Sn.j  +  t(Xr) 

Rectangular 

Sr  -  *„_!+!  ( 3X„  -  Xr.!) 

Implicit  Adams 

Second-Order 

Sr  -  Xr.j  +  t(2Xr  -  Xr_i) 

— 

SECTION  VII 


A  T/2,  T/3  CLOSED-LOOP  SIMULATION  CASE  STUDY 

A.  INTRODUCTION 

The  topic  of  this  section  is  a  case  study  using  the  tools  developed 
in  Sections  III  through  VI.  The  rates  involved  in  this  closed-loop 
analysis  force  the  use  of  switch  decomposition.  Later  sections  will 
treat  case  studies  where  switch  decomposition  can  be  avoided. 

Specifically,  we  investigate  a  case  where  an  idealized  control  ele- 
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ment  (1/s  )  is  under  the  influence  of  a  continuous  feedback  controller 
with  an  (idealized)  compensation  network  in  the  forward  path.  The 
objective  is  to  generate  Bode  plots  given  (1)  a  digital  implementation 
of  the  continuous  controller  and  (2)  a  digital  simulation  of  the  overall 
system;  treating  both  the  continuous  and  discrete  controller  cases.  Of 
special  interest  is  the  case  where  the  compensation  is  modeled  on  one 
computer  and  the  controlled  element  is  modeled  on  a  different  computer, 
with  each  computer  working  in  a  different  frame  time.  In  all,  four 
cases  will  be  discussed. 

B.  PROBLEM  DESCRIPTION 

The  situation  described  by  Case  I  in  Figure  40  assumes  a  controlled 
2 

element  (1/s  )  under  the  influence  of  a  continuous  feedback  controller 
with  an  (idealized)  compensation  network  2(s  +  1)  in  the  forward  path. 
Case  II  depicts  the  same  controlled  element  under  the  influence  of  a 
discrete  feedback  controller  which  smooths  the  output  of  the  digital 
computer  with  a  ZOH  (passing  on  a  "staircased"  signal  to  the  control 
point).  The  discretized  version  of  2(s  +  1)  was  computed,  using  the 
first  back  difference  algorithm  resulting  in  the  model  (42z  -  40)/z  (at 
a  sampling  rate  of  20  Hz)  . 

2 

Case  III  is  a  simulation  of  Case  II  (or  I)  wherein  the  plant,  1/s  , 
is  modeled  using  a  "substituticn-for-s"  rule  (Tustln  transform) .  The 
output  variable  C  is  modeled  as  the  output  of  a  storage  register  (ZOH) . 


Finally,  Case  IV  depicts  a  situation  wherein  one  part  of  a  simula¬ 
tion  is  coded  for  one  computer  while  another  part  is  coded  on  a  second 
computer.  Typically,  the  computers  are  working  in  different  frame  times 


and  therefore  will,  on  occasion,  pass  "old  data"  back  and  forth.  It  is 
assumed  that  the  compensation  is  modeled  in  a  0.05  sec  time  frame,  while 
the  plant  is  modeled  in  a  0.075  sec  time  frame.  Data  transfer  between 
the  two  computers  is  via  appropriate  buffer  registers,  modeled  as  a  ZOH 
in  a  T/3  time  frame  (M^)  and  a  ZOH  in  a  T/2  time  frame  (M2). 

This  completes  the  problem  description;  we  may  now  write  the  appro¬ 
priate  equations  for  each  case  and  discuss  the  analytical  difficulties. 

Case  I  is  straightforward  since 

C(s)  -  [i  +  G1G2]-1  G1G2R  (214) 

Case  II  appears  straightforward  but  does  present  a  contradiction  in  the 
output  equation 


CT  *  [i  +  (GjM)T  CT]  (g^m)T  GTRT  (215) 


That  is. 


(GiM)T 


ll  (z  +  *> 

2  (z  -  l)2 


(216) 


is  readily  computed  using  a  transform  table.  However,  it  was  not  the 
intent  of  the  designer,  who  used  a  substitution  technique  to  model 
2(s  +  1)  as  (42z  -  40)/z,  to  compute 

G*  -  [2(s  +  1)]T  (217) 

T 

using  the  z-transform.  Indeed,  what  does  [s]  mean?  Clearly,  the 

T 

intent  of  the  simulation  designer  was  to  assign  G2  *  (42z  -  40) /z. 

> 


This  difficulty  surfaces  again  in  Case  III  since 


CT  .  [l  +  gJgT]-1  gTGTrT  (218) 

T 

after  taking  due  note  that  M  =  1. 

T  T 

Now  it  is  necessary  to  interpret  both  G^  and  G^  as  given  functions 
of  z  rather  than  z-transform  operations.  For  example,  assign  (via  the 
Tustin  transform) 


=1 


2  2 

if  U-±.  n 

4  (z  -  l)2 


(219) 


rather  than  (via  the  z-transform) 


Tz 

(z  -  l)2 


(220) 


Difficulties  in  the  "assignment"  procedure  also 
write  the  output  equation  for  Case  IV  (see  Figure  41). 


surface  when  we 


Figure  41.  Case  IV,  Switch  Decomposition  Model 
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In  a  similar  manner 


(W2*M3W3)t 


T 


l[_L 

il  - 

1  j  e-sT/3 

e-s2T/3|  j 

(  .  e«T/2. 

s 

T 

/l  -  e- 

(T/3 

1 

e-aT/3 

e“s2T/3l\ 

\  ’ 

esT/2 

esT/6 

e-st/6  \J 

(225) 


Clearing  through  and  writing  Equation  225  in  terms  of  the  advanced 
z-transform  for  1/s: 


(w2*m3w3>T 


1  e-3Tea2T/3 

e-sTes 2T/3  e-sTesT/3 

e-8TgST/3 

ro 

1 

CO 

H 

s  s 

s  ""  s 

s 

S 

qsT / 2  gsT/6 

esT/6  e-sTes5T/6 

esTes5f / 6 

e-sTesT/2 

s  ~  s 

s  s 

s 

s 

1  0  0 

0  1  0 


(226) 


However,  the  meaning  of  (W3*G2W3)T  or  (W2*G3W2)T  Is  not  clear.  When  the 
designer  specified  G2  ”  (42z  -  40) /z,  he  intended  (for  a  T  time  frame) 
the  recursion  relationship 


In  a  like  Banner,  the  Tustln  substitution  for  1/s  was  meant  to  yield 


X4(n)  -  2X4(n  -  1)  -  X4(n  -  2)  +  j?  [X3(n)  +  2X3(n  -  1)  +  X3(n  -  2)] 


(228) 


Clearly,  the  Information  available  In  Case  IV  Is 


GT/3  .  ‘ill..- A  0 

2  z 


,  T/3  -  0.05  ,  z  «  esT/3 


(229) 


16  (z  -  1)' 


,  T/2  -  0.075  ,  z  =  esT/2 


(230) 


T 

Using  only  the  given  computer  code,  how  does  one  compute  (W3  G2W3^  and 

■»  T  ^ 

(W2^giW2^  -  This  problem  Is  treated  next. 

C.  DECOMPOS IT IOH  OF  G*/Z  AHD  G*/3 


Using  the  results  of  Section  VI,  write 


(w3*g2w3)t  -  (w3ykcl/3w3) 


(231) 


This  Is  a  simple  but  crucial  step  because  we  do  not  know  (or  care)  what 

T/3 

the  underlying  G(s)  was.  However,  we  are  given  Gj  as  computer  code 

or,  equivalently,  a  pulse  transfer  function  in  a  T/3  time  frame.  The 
T/3 

G3  transfer  function  from  the  Case  IV  simulation  example  has  the  form: 


a  z  +  a, 
o  _ l 


-  esT/3 


(232) 


Equation  231  has  the  terms 


T /3  T 

(w3*G2/  W3) 


z-3[22G5/3] 


z-3[2G|/3] 

z~3[z2G2/3 


(233) 


where  z  ■  esT/3.  The  final  result  will  be  in  a  T  time  frame  and  there¬ 
fore  the  z-3  in  the  superdiagonal  terms  can  be  factored  out  and  treated 
as  a  z~l,  where  z  ”  es^. 

Only  the  entries  for  the  first  column  need  be  computed: 


aQz  +  a^ 

T 

z(a0z  +  ai) 

T 

z2(a0z  +  31) 

z 

9 

z 

9 

z 

(234) 


Consider  the  first  term: 


a0z  +  ax 

T 

=  j 

faoP  +  al  z  dp 

z 

2n  j  J 

P  z  -  p3  P 

_d 

(aoP  +  al)z 

dp 

Z  -  p3 

p=0 

(235) 


(236) 


This  result  is  easily  verified  using  the  algorithm  of  Section  IV. 
In  the  T/3  time  frame,  algebraic  long  division  gives 


aQ2  +  al 


aQ  +  aiz“l  +  Oz~^  +  0z~3  + 


(237) 


The  root  at  z  -  0  in  the  T/3  time  frame  remains  at  z  -  0  in  the  T  time 
frame.  Therefore, 


a0z  +  ai 


ciz  +  c0 


(238) 


where 


ao  0 


0  0 


(239) 


and  we  see 


a0z  +  ai 


(240) 


The  second  term  in  Equation  234  can,  of  course,  be  evaluated  using  resi¬ 
due  theory: 


([•o*  +  ali) 


1  f  (aoP  + 

2l,iJT  p(p  - 


ai)z 

p3)  dp 


(aPP  +  ai)z 


z  -  p3 


p“0 


al 


(241) 


and 


r  j  /  (aQp  +  ai)z 

[z2(a0z  +  a^JJT  -  /  - r —  dp  i  0  (no  poles  enclosed) 

z  -  pJ 


(242) 


In  this  regard,  note  that  the  algorithm  of  Section  IV  was  set  up  to 
treat  only  rational  functions  since  polynomials  in  z  can  be  easily 
treated  by  inspection.  That  is: 
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The  first  tern  of  column  one  Is 


T2(z  +  l)2 

T  1  f  T2  (p+1)2  _ 

z  dp 

16(z  -  l)2 

2*J  JT  16  (p-l)2  z_p2  p 

res 

.  I2  ,  T2  (£±1)2z  , 

-  i2 

1  +  ~ 

16  +  16  dp  p(z-p2) 

16 

P-1 

( 

T2  (z2  +  6z  +_1) 

16  (z  -  l)2 

»  —  — 

The  second  tern  of  the  first  colunn  Is  easier  to  compute: 


1  f t£  p(p  +  1)2  2  dp 

2*j  J  16  (p  -  l)2  (z  -  p2)  P 


res 


P-1 


T2 

16 


St¬ 

rip 


z(p  +  l)2 
(z  -  p2) 


P-1 


.  z(z  +  i) 

4  (z  -  l)2 


Again  use  the  push  down  column  principle  and  write 


0»2*<W 


fg  (z2  +  62  +  1) 

j2 

*7  z(z  +  1) 

A 


t2 

~  (z  +  1) 

<z2  +  6z  +  1) 


(z  -  l)‘ 
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The  algorithm  of  Section  V  can  also  be  used.  To  demonstrate,  derive  the 
first  entry  by  observing 


t£  (z  ±  l)2 
16  (z  -  l)2 


[l  +  4z_1  +  8z~2  +  12z-3  +  16z-4  +  •••]  (252) 


Since  a  root  at  z  -  1  in  T/2  remains  at  z  -  1  in  a  T  time  frame,  write 


T 2  fz  4-  1)2 

16  (z  +  l)2 


C2Z  +  C[z  +  c0 
z2  -  2z  +  1 


,sT 


(253) 


the  numerator  coefficients  are  easily  found 


c2 

T2 

ci 

_C°. 

“  16 

L 

1 

8 

16 


0 

1 

8 


1 

-2 

1 


16 


(254) 


checking  the  residue  results. 
D.  RESULTS 


With  the  assignments  of  the  previous  subsection  all  terms  in  Equa¬ 
tion  221  are  defined  and  the  Bode  plots  descriptive  of  Cases  I-IV  can  be 
computed.  The  results  (magnitude  plots  only)  are  shown  in  Figures  42 
through  47. 

Inspection  of  Figure  42  discloses  no  surprises.  The  digitally  con¬ 
trolled  system  is  a  reasonably  faithful  reproduction  of  the  analog 
system  until  the  folding  frequency  (approximately  equal  to  62.8  rad/sec) 
is  passed.  Notice  that  in  the  discretely  controlled  system,  minimum 
response  points  In  the  Bode  plot  (notches)  occur  at  multiples  of  sampl¬ 
ing  frequency  (approximately  equal  to  125.66  rad/sec) . 

The  comparison  of  the  single-rate  simulation  against  the  baseline 
design  exhibits  fidelity  over  a  shorter  low-frequency  range  (Figure  43) . 


V 


Figure  43.  Magnitude  Plot;  Case  I,  Case  III 


Of  particular  importance  is  that  the  aliased  bands  exhibit  a  much  higher 
amplitude  response  than  did  the  aliased  bands  of  the  digitally  con¬ 
trolled  system  (Case  II).  Moreover,  the  "notches"  now  occur  at  multi¬ 
ples  of  the  folding  frequency  rather  than  the  sampling  frequency.  In 
addition,  there  are  very  sharp  notches  which  occur  close  to  odd  multi¬ 
ples  of  the  folding  frequency;  these  are  a  consequence  of  the  zeros  of 
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the  Tustln  transform  introduced  by  (z  +  1)  .  A  direct  comparison 
between  Cases  II  and  III  is  given  in  Figure  44. 

Figure  45  compares  the  two  independent  processor  case  (Case  IV) 
against  the  continuous  baseline  design.  Large,  sharp  resonant  peaks 
have  been  introduced  in  the  aliased  bands  and,  in  addition,  there  is  a 
significant  overshoot  in  the  first  fold.  A  comparison  between  the  two 
rate  simulation  (Case  IV)  and  the  digitally  controlled  continuous  system 
(Case  II)  is  given  in  Figure  46.  Figure  47  compares  all  the  cases. 

There  are  significant  differences  in  the  spectral  content  of  the 
four  cases  which  would  be  hidden  if  one  only  looked  at  the  sampled  spec¬ 
trum  (that  is,  looked  only  at  the  frequency  content  from  zero  to  the 
first  folding  frequency).  Even  in  the  first  fold,  there  is  a  signifi¬ 
cant  difference  in  the  Bode  plot  of  the  continuous  case  and  the  two-rate 
simulation;  the  reason  for  the  added  overshoot  in  the  two-rate  simula¬ 
tion  will  be  discussed  in  the  next  subsection. 

E.  INTRODUCTION  OF  LIGHTLY  DAMPED  MODES 

We  may  use  the  illustrative  example  of  Section  V  to  discuss  the 
additional  modes  introduced  by  a  multi-rate  process.  Specifically,  for 
Figure  48,  let 
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esT/6 


M  -  3  ,  N  -  2 
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Case  32  • 


Figure  48.  A  Phantom  Sampler  Formulation  of  a 
T/N,  T/M  Sampling  Format 


so,  as  found  in  Section  V, 


cT/3 


_ 22fz2+e~5T/61 _ 

(z_e-2T/3j(z_e-T/3)[z2+  ( e-T/3 )z+e-2T/3 j 


2  =  esT/3  (258) 


Observe  that  the  T/2,  T/3  sampling  format  has  produced  additional 
lightly  damped  modes  in  the  output  response.  The  reality  of  the  addi¬ 
tional  modes  can  be  better  appreciated  by  first  plotting  the  continuous 
variable  C(t)  and  then  picking  off  the  T/3  sample  points.  This  is  done 
in  Figure  49.  Joining  the  sample  points  with  a  smooth  curve  emphasizes 
the  lightly  damped  nature  of  the  response.  This  effect  was  also  present 
in  the  two-rate  simulation  analysis  of  the  previous  subsection  (recall 
the  additional  "overshoot"  in  the  first  fold  of  Figure  45)  • 

F.  SECTION  SUMMARY 

A  case  study,  which  required  the  use  of  switch  decomposition,  was 
used  to  demonstrate  the  significant  spectral  differences  that  occur  when 
a  closed-loop  system  (either  an  analog  or  digital  controller)  is  com¬ 
pared  with  an  all-digital  simulation  of  the  closed-loop  system.  The 
example  also  treated  the  problems  encountered  when  a  simulation  software 
is  put  up  on  two  different  (Independent)  computers,  each  working  in  a 
slightly  different  frame  time. 

In  following  sections,  simulation  case  studies  amenable  to  more  dir¬ 
ect  analysis  (e.g.,  no  need  for  switch  decomposition)  will  be  studied. 


Ill 


SECTION  VIII 


INTERPRETATION  OF  MULTI-RATE  FREQUENCY  RESPONSE 


A.  INTRODUCTION 

A  simulation  case  study  of  the  previous  section  required  the  use  of 
switch  decomposition.  This  section  treats  an  open-loop  case  which  can 
be  analyzed  by  the  more  direct  "scalar"  approach  discussed  in  Section  V, 
affording  the  opportunity  to  gain  insight  in  interpreting  the  frequency 
response  of  a  multi-rate  system. 

It  is  shown,  for  a  T/2  input,  T/3  output,  that  three  components  are 
required  to  produce  an  exact  steady-state  match  to  the  output  sample 
sequence . 

This  simple  case  study  also  affords  the  opportunity  to  contrast  the 
mathematics  of  the  scalar  approach  against  that  of  vector  switch  decom¬ 
position.  This  is  done  by  working  the  scalar  example  through  using 
switch  decomposition  notation. 

B.  A  T/2  INPUT,  T/3  OUTPUT  SIMULATION  EXAMPLE 

The  study  is  defined  in  Figure  50.  Our  objective  to  verify  that 

three  sinusoidal  components  are  needed  to  exactly  match  the  steady  state 
T/3 

sampled  points,  C  .  From  Figure  50,  with  r  -  sin  bt,  write 


CT/3  -  |m2 


1  -  e-T/2 
z  -  e“T/2 


z  sin  bT/2 


z^  -  2(co8  bT/2)z  +  1 


|T/3 


»sT/2 


(259) 


A  term  of 


(260) 


2 


2 


1/5 


Figure  50.  A  Two  Rate  Experiment 


is  in  a  T/2  time  frame.  Since  the  commensurate  time  frame  is  T/6,  z' 

T  /  6 

must  be  substituted  for  z.  Moreover,  M2  becomes,  in  a  T/6  time  frame 


T/6 

1  -  g-sT/2 
s 


23  -  1  Z  _  2?  +  z  +  1 

z3  z  -  1  =  z2 


(261) 


Therefore 


-T/3  = 


+  z  +  l 


( z 


{ 1  ,  e-T/2)  z3  sln  bT/2 
3  -  e~T/2j(z6  _  2  cos  z3  +  l) 


T/3 


(l  -  e^T/2]  sin  ^  (z^  +  +  z) 

z9  .  (e-T/2  +  2  cos  ^§)z6  +  (l  +  2e-T/2Cos  ^|j  x3  -  e"T/2 

(262) 

As  an  alternative  to  the  use  of  Section  IV  for  the  evaluation  of  Equa¬ 
tion  262,  we  will  use  Sklansky's  Identity  (Appendix  B)  and  a  frame  time 
of  T  *  0.1  seconds: 


[GT/2(z)]T  -  [G(z)  +  G (-z ) ] 


(263) 
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Since  the  ratio  of  T/6  to  T/3  is 
given  by 


2,  the  rational  form  of  C^/3  be 


cT/3 


Rl6Z^  +  Risz^  +  •••  +  Rnz  +  Rio 

R29Z9  +  R26z6  +  r23z3  +  r20 


(264) 


where  the  R^j  coefficients  are  as  listed  in  Table  3. 

Dividing  the  denominator  of  Equation  264  into  the  numerator  gives 
the  transient  response  in  a  T/3  time  frame.  This  response  is  sketched 
in  Figure  51  as  a  solid  line  but  the  actual  discrete  response  follows  a 
2— 1—2— 1—2— 1—2 . . .  repetition  pattern  (note  the  insert  in  Figure  51). 

Returning  to  Equation  264,  observe  z2  -  (2  cos  bT)z  +1,  in  a  T  time 
frame,  is  a  factor  of  the  denominator.  After  factoring  out  this  term, 
the  denominator  has  the  form: 


(R29z9  +  R26z6  +  R23z3  +  r2o) |z=esT/3 

=  R29z3  +  R26Z2  +  R23z  +  r20 | z_esT 


=  [z2  -  (2  cos  bT)z  +  1 ] (z  -  .904837421)  (265) 


Write  Equation  264,  recognizing  Riq  =  0,  as 


cT/3 


(Rj6Z6  +  R15Z5  +  »•»  +  Rn) 
[sin  (bT)]z2 


z»e 


sT/2  z  ~  .904837421  |z_esT 


z  sin  bT 


.2  _ 


(2  cos  bT)z  +  1 


115 


- --  --  Cf  — » 


(266) 


Figure  51 


Versus  Time 


Equation  266  has  the  Input  factored  out,  in  a  T  time  frame,  and 
therefore  is  equivalent  to  the  switch  decomposition  approach. 

The  spectral  coefficients  are  determined  from  the  equation 


A  +  IB 

n  *■'  n 


_L 


+  R \ ^2^  +  ♦♦  +  R 


U 


n  n  3  sin  bT 


z-Uo)bT/2 


(267) 


z  -  .904837421  li4u>nT/2 


where 


b  + 

T 


n  -  0,  1,  2 


(268) 


and 


[CT/3] 


ss 


2  (An  sln  “n*-  +  Bn  cos  “n*-) 
_n«0 


T/3 


(269) 


Notice  an  important  point  —  the  scale  factor  sin  bT  and  the  coeffi¬ 
cients  R j  j  through  R^  are  held  fixed  at  the  input  frequency  but  the 

various  z's,  in  their  respective  frame  times,  are  evaluated  at  wn!  A 

consequence  of  ignoring  this  point,  and  running  the  program  from  scratch 

each  time  (generating  R^  *  Ri6>  etc.)  in  terms  of  u>n  Instead  of  b,  will 
result  in  an  erroneous  set  of  coefficients.  To  verify  this,  we  imple¬ 
ment  both  procedures  and  compare  the  results. 

In  Table  4,  the  coefficients  (and  sin  bT)  are  first  evaluated  using 
the  input  frequency  b.  They  are  then  "frozen"  and  the  next  two  spectral 
components  are  generated.  In  Table  4b,  the  evaluation  is  repeated, 
n  ■  0,  1,  2,  treating  as  a  new  input  frequency.  As  can  be  seen,  the 
first  and  third  coefficients  are  unchanged  but  the  second  components 
disagree. 


We  next  generate  the  sum  of  the  three  sine  waves  for  each  case  and 
see  how  they  compare  against  the  actual  sinusoidal  transient  response. 


TABLE  4.  SPECTRAL  COEFFICIENTS  FOR  EQUATION  277 


MAG  -3.010199322 

-3.010199412  MAG 

n-0  -47.39926233 

-47.39926170  4. 

u>  1.000000000 

1.000000000  w 

-43. 30234384 

-35.68901104 

n-1  -287.3992625 

-242.4231347 

63.83185308 

63.83185308 

-43.38592638 

-43-38592638 

n-2  -707.3992644 

-707.3992644 

126.6637062 

126.6637062 

(a) 

(b) 

as  tabulated  in  Table  5.  It  can  be  seen  that  the  Table  4a  coefficients 
generate  the  sums  of  three  sine  waves  which  agree  extremely  well  with 
the  transient  response  at  the  T/3  sampling  instant.  The  Table  4b  set 
produces  a  sum  of  three  sine  waves,  which  is  far  less  precise  (sometimes 
disagreeing  in  the  first  significant  figure). 

The  important  result  is  that  the  T/3  output  signal  required  three 

sinusoidal  components  for  the  correct  steady-state  fit.  Contrasting 

this  result  with  the  switch  decoraposition/Bode  approach  of  Section  V 

affords  an  interesting  comparison.  First  of  all,  the  approach  of  Sec- 

T 

tion  V  is  "dimensioned"  by  (W *R)  ,  which  in  this  case  is  two  since  R  is 
sampled  in  a  T/2  time  frame.  How  is  this  reconciled  with  the  present 
result,  which  insists  on  three  output  components? 

In  reality,  there  is  no  disagreement,  since  the  Section  V  approach 
would  use  two  Bode  plots  to  pick  off  the  same  three  components.  We  can 
verify  this  by  reworking  the  subject  example  using  switch  decomposition. 


TABLE  5.  COMPARISON  OF  STEADY-STATE  TRANSIENT  RESPONSES 


—  69393 

-0.69394 

-0.68198 

-.69393 

-0.69394 

-0. 70351 

-.68625 

-0.68626 

-0.68810 

-.67686 

-0.67687 

-0.66536 

-.67686 

-0.67687 

-0.68721 

a) 

Transient  response  of^g 

Equation  264,  from  z„ 

.  -397  ’  sT/3 

to  z  ,  where  z**e  , 

-.66578 

-.65304 

-0.66579 

-0.65305 

-0.66638 

-0.64209 

-.65304 

-0.65305 

-0.66404 

T  -  0.1. 

-.63866 

-0.63867 

-0.63800 

-.62268 

-0.62269 

-0.61240 

b) 

Sum  of  three  sine  waves 

-.62268 

-0.62269 

-0.63424 

using  amplitudes  and 

-.60515 

-0.60516 

-0.60325 

phase  angles  of  Table  4a 

-.58611 

-0.58612 

-0.57660 

-.58611 

-0.58612 

-0.59810 

c) 

Sum  of  three  sine  waves 

-.56560 

-0.56561 

-0.56247 

using  amplitudes  and 

-.54368 

-0.54369 

-0.53503 

phase  angles  of  Table  4b 

-.54368 

-0.54369 

-0.55598 

-.52040 

-0.52041 

-0.51606 

—  49582 

-0.49582 

-0.48812 

—  49582 

-0.49582 

-0.50831 

-.47000 

-0.47000 

-0.46451 

—  44300 

-0.44301 

-0.43633 

-.44300 

-0.44301 

-0.45556 

-.41490 

-0.41490 

-0.40831 

-.38576 

-0.38576 

-0.38018 

-.38576 

-0.38576 

-0.39826 

—  35566 

-0.35566 

-0.34803 

-.32466 

-0.32466 

-0.32023 

-.32466 

-0.32466 

-0.33698 

—  29286 

-0.29286 

-0.28427 

-.26032 

-0.26032 

-0.25708 

-.26032 

-0.26032 

-0.27233 

—  22714 

-0.22713 

-0.21767 

—  19338 

-0. 19338 

-0.19137 

—  19338 

-0. 19338 

-0.20496 

-.15914 

-0. 15914 

-0.14890 

—  12451 

-0.12450 

-0.12374 

—  12451 

-0.12450 

-0.13555 

C.  SWITCH  DECOMPOSITION  DEVELOPMENT 


For  comparison  purposes ,  we  next  use  switch  decomposition  on  the 
subject  example  and  redraw  Figure  50  as  Figure  52.  From  Figure  52, 


CT/3  -  W3(w3jkM2W2)T(w2jk  s-^-T  M2W2)T(W2*R)T 


(270) 


» 

o 

_ 1 

e-T/2_e-T  i  _  e-T/2 

„  [i,  e-sT/3,  e-(2/3)sT] 

1  0 

zfl-e-T/2)  e-T/2_e-T 

I 

H 

o 

_ 1 

z  -  e-^- 

RT 

(esT/2R 


The  theory  of  Section  V  gives  the  coefficients  as 


(271) 


An  +  JBn 


l  +  e-sT/3  j  e-(2/3)sT ] 
3 


x 


'  e-T/2  -  e-T  }  _  e-T/2  “ 

i 

(l-e-T/2),  e~T/2  _  e-T 

l$.bT/2 

-T 

z  -  e  A 

(272) 


n  -  0,  1,  2 


Figure  52.  Switch  Decomposition  Model  for  Figure  48 


Therefore 


U) 


n 


b 


+ 


2nn 

T 


[l  +  li-^T/3  !  1  4-~(  2/3)to  T  ] 
An  +  JBn  -  - T - 


(273) 


e-T/2  _  e-T  1  _  e“T/2 


L(l-e"T/2)4n,I  e-T/2  -  e-T 

1 4%T - e^ 


1 

UbT/2 


n  =  0,  1,  2  (or  n  =  0,  ±1) 

From  Equation  273,  it  is  evident  that  there  is  no  conflict  between 
switch  decomposition  and  the  scalar  approach  since,  for  n  *  0,1,2,  we 
obtain  the  spectral  coefficients  listed  in  Table  4a.  The  switch  decom¬ 
position  approach  confirms  that  three  components  are  required  (n  =  0, 
±0.1);  in  addition,  (W*R)T  requires  that  the  input  remain  fixed  at 
b  rad/sec.  The  tradeoff  between  the  scalar  approach  and  switch  decompo¬ 
sition  is  now  clear.  Switch  decomposition  has  a  greater  dimensionality 
problem,  but  it  also  has  a  format  which  protects  the  user  from  making  an 
error  since  the  correct  uses  of  u)n  and  b  are  explicitly  called  out*  In 
this  regard,  the  user  must  have  a  clear  understanding  of  the  scalar 
approach  in  order  to  use  u>n  and  b  in  the  correct  sequence. 

D.  SKCTIOD  SUMMARY 

A  scalar  example  was  used  to  compare  the  "scalar"  approach  with 
switch  decomposition.  It  was  shown  that  both  yield  identical  results. 
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SECTION  IX 


A  THREE-BATE  SIMULATOR  CASE  STUDY 


A.  INTRODUCTION 

The  simulation  error  analysis  case  studies,  up  to  this  point,  have 
been  restricted  to  low-order  examples.  Here  we  depart  from  this  motif 
by  investigating  a  simulation  involving  three  computers,  each  working 
within  a  different  frame  time. 

As  we  shall  see,  no  new  theoretical  tools,  other  than  those  already 
introduced,  will  be  needed  to  carry  out  the  analysis.  However,  the 
reader  is  alerted  to  the  fact  that  many  computational  difficulties  were 
encountered  during  the  course  of  this  case  study.  These  difficulties 
forced,  for  example,  the  "invention"  of  the  multi-rate  algorithm  given 
in  Section  IV.  The  nature  of  the  difficulties  encountered  will  be  indi¬ 
cated  by  first  describing  a  lower  dimensional  example  which  retains  the 
basic  structure  of  the  large-scale  system  study.  After  this,  we  proceed 
with  the  analysis  of  the  title  study. 

B.  A  SIMPLIFIED  THREE-RATE  STUDY  AND  ASSOCIATED 

DIMENSIONAL  DIFFICULTIES 


For  clarity,  utilize  a  precise  set  of  s-plane  parameters  in  order  to 
produce  numerically  convenient  numbers  in  the  z-plane.  Consider  the 
continuous  system  of  Figure  53a  which  is  to  be  simulated  with  the  three 
rate  configuration  of  Figure  53b).  In  Figure  53,  let  T  ■  0.1  sec 
(10  Hz), 


4.52569504 
s  +  4.52569504 


2.2222222 

s  +  2.2222222 


(274) 


Figure  53.  A  Three  Rate  Example 


and 


M 


1 


t  _  e-aT/16 
s 


m2 


(275) 


Assume  the  "digitization"  required  in  Figure  53b  is  carried  out 
using  rectangular  integntion  [s  +  (z-l)/Tz].  For  comparison  purposes, 
the  transfer  function  of  the  continuous,  baseline  system  is 


£ 

R 


20^2 

1  +  GjG2 


_ 20.  11420018 _ 

s2  +  6.747917262s  +  20.11420018 


(276) 
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Utilizing  the  analysis  technique  of  Section  III  gives  C 


CT/3  .  [  2M2G2^  [l  +  (Mic{^16M2^I6)T^2G21^2]  (MiGi216RT^6) 


[«20l/2(HlGl/16R^‘6)  ] 


Step  I:  Begin  the  evaluation  of  Equation  277  with  the  term 


(HlGyl6M|/16)T/2  -  (hI/I6gI/16«I/16)T/2 


(278) 


Therefore 


(GT/16MT/16}T/2 


z  -  .972492473 


,8  z  -  1 


since  rectangular  integration  gives 


i  aT 

1  +  aT^z 

s  -  ( - — ) 

M  +  aT' 


-sT/2 ' 


Rearrange  Equation  279  as 


[gI/16m!/16] 


z  -  1  -8 


sST/16 


(280) 


(281) 


k0(z7  +  z6  +  z5  +  +  l) 

z  7  -  k  ^z® 


(282) 


Using  the  techniques  of  Section  IV,  divide  the  denominator  into  the 
numerator  and  save  every  eighth  number,  since  the  ratio  between  the  two 

g 

frame  times  is  eight  (see  Table  6,  and  also  note  0®  -  0;  kj  -  0.8). 


test;--,  ‘it- * 

. .  c  '  * 


TABLE  6.  THE  "a"  VECTOR,  STEP  I 


a7 

0.027507527 

a6 

0. 194498493 

a5 

0.155598794 

a4 

0.124479035 

a3 

0.099583228 

a2 

0.079666583 

al 

0.063733266 

ao 

0.050986613 

Setting 


[G^16  Mj/16 


T/2 

] 


C72^  +  052^ 

(z  -  . 


+  ♦  ♦  • 

8)z6 


+  CQ 


g 


Solving  the  set  of  Equation  284  gives  the  T/2  frame  time  transfer 
function 


[„T/16.,T/16 1T^2  .027507527Z  +  .172492471 

l°l  "2  1  - (Tr~. 8) - 


RqZ  +  Rj 
z  -  •  8 


z.esT/2 


(285) 


In  Equation  284,  only  the  first  two  columns  of  the  A  matrix  are 
needed  since  the  "C"  vector  has  only  two  components.  Further,  note 


c5  “  a5  ~  *8a6  =  0 


c 


4 


=  a4  ~  •  8a^  ■=  0 


c0 


0 


(286) 


verifying  that  z  is  a  common  factor  of  both  the  numerator  and  denomi¬ 
nator  of  Equation  283. 

At  this  point  it  is  clear  that  modeling  the  "buffer  registers" 
between  computers  as  ZOHs  introduces  dimensionality  problems  which  get 
larger  as  the  ratio  between  the  frame  times  increases.  The  evaluation 
of  Equation  283,  using  residue  theory,  would  be  tedious  since  it 
involves  a  root  at  the  origin  of  multiplicity  six!  This  multiplicity 
would  increase  significantly  if  the  rates  had  been,  for  example,  20  Hz 
and  156  Hz  since  the  ratios  between  frame  times,  in  the  last  step,  jumps 
to  39  (see  examples  A-3  through  A-5  of  Appendix  A).  While  the  computa¬ 
tional  burden  would  become  simply  enormous  using  residue  theory,  the 
algorithm  of  Section  IV  is  affected  in  a  far  less  sensitive  way;  since 
the  only  change  required  is  the  storage  of  every  39th  value  (instead  of 
every  eighth  value)  of  the  power  series  expansion. 
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A 


Step  II: 


Since 


gT/2  2.22222  '  .lz 

2  '  Is  +2-2222\,..i-l)/(T/2U  2  -  '9 

the  evaluation  of 

[•j-1  -  i  +  (m1gT/16m5/16)  G^2  =  1  4 


»  T  -  .1  (287) 


• 1z(RqZ  +  Rj ) 

(z  -  .8)  (z  -  .9) 

(288) 


_ (z  -  . 8) (z  -  .9) _ 

(1  +  . 1R0) z2  +  (.1RX  -  1. 7)z  +  .72 


(289) 


is  routine. 


Dividing  out  the  lead  coefficient  of  the  denominator: 


[•r1 


. 99 7256793 (z-.8)(z-. 9) 
zz  -  1.678134619z  +  .718024891 


R^z-. 8) (z-. 9) 
z^  +  R^z  +  R^ 


(290) 


Computing  2G^  is  also  straightforward: 


2g|/2  f • J _1 


■  2R,(z  -  .8)z 


z^  +  R^z  +  Rj 


(291) 


Step  III: 

Next,  treat  the  input  term,  assuming  a  unit  amplitude  sine  wave  at 
b  rad/sec: 


-  [g(/,6RT/16]1 


since  is  unity 


(292) 


. 027507527z 
z  -  .972492473  *  JT 


z  sin  bT/16  ]T/2 

2  cos (bT/16)z  +  -J;_eST/16 


In  the  computer  program,  b  is  left  as  a  free  input  parameter.  Here  the 
interest  is  in  scoping  dimensionality  problems  associated  with  the  eval¬ 
uation  of  Equation  277;  therefore  we  pick  a  numerically  convenient  value 
for  b  in  order  to  obtain  an  explicit  expression  for  Equation  292. 

Let 


bT 

COS  —T 


(293) 


giving 


[®T 


T/16„T/16 


T/2 

] 


(. 02750752 7) ( . 08035Q925)z2 


(z  -  .972492473) (z2  -  1.993533274z  +  1) 

,T/2 


T/2 


Z3  +  RyZ^  +  RgZ  +  Rg 


In  the  T/2  time  frame,  the  poles  of  Equation  294  map  into 


(294) 


(z  -  .8)  (z*  -  1- 6z  +  1) 


z3  -  2.4z2  +  2. 28z  -  .8 


(295) 


giving  the  form  for  a 


T/2 


[gI/16rt/16) 


T/2  3  2 

'  C3Z  +  C2Z  +  cjz  +  eg 

z3  -  2.4z2  +  2. 28z  -  .8 


(296) 


Dividing  the  denominator  into  the  numerator  of  Equation  294,  gives  the 
"a"  vector  (save  each  eighth  number).  Therefore 


■  V 

a3 

0 

0 

0  " 

‘  1  ' 

c2 

a2 

a3 

0 

0 

-2.4 

C1 

al 

a2 

a3 

0 

2.28 

-C0. 

.  ao 

al 

a2 

a3. 

-.8 

(297) 


Solving  Equation  297  gives  the  solution;  see  Table  7  for  a 
of  the  "a"  and  "c"  vectors.  Explicitly, 

rGT/16RT/16iT/2  m  (.071804007Z  +  .Q5238627)z 
1  (z  -  .8)(z2  -  1.6z  +  1) 

(^10z  +  ^ll^z 

(z  -  .8) (z2  -  1.6z  +  1) 


TABLE  7.  THE  "a"  AND  "c"  VECTORS,  STEP  III 


a3 

0.000000000 

C3 

0.000000000 

a2 

0.071804007 

C2 

0.071804007 

al 

0.224715895 

C1 

0.052386278 

ao 

0.375604980 

co 

0.000000000 

Step  IV: 

The  last  step  introduces  another  dimensionality  problem, 
result,  to  this  point,  is 


CT/3  . 


MacriMiorV"”) 


T/16  T/16 


T/2 


T/3 


listing 


(298) 


The 


,-sT/2  ( .  2R3)  j[z- — "TSjz 


(Ripz  +  Rll)z 


T/3 


1 

/(l  -  e-sT)\T/6 _ (r12z  +  R13)z2 

\  a  +  Rj^z^  +  +  Rj^z  +  ^17 

J9 


z3  -  1  z  r12z3  +  r13 

"  0  —  ■■ —  •  ■ ■  - 

z^  z  -  1  z12  +  Rj^j9  +  R^z^  +  Rigz^  +  Rj^ 


(300) 


(301) 


Cancelling  powers  of  z  and  dividing  z  -  1  into  z  -  1  gives 


(z^  +  z  +  i)z^(r^2z3  +  R l 3 ) 

*12  +  Rj^z^  +  Rjjz^  +  R^z^  +  Rj^ 


(302) 


Considering  that  one  starts  with  two  first-order  transfer  functions 
in  the  s-plane,  the  ninth/twelfth  order  transfer  function  in  z  is  dis¬ 
concerting.  It  is  a  good  example  of  how  the  order  of  the  system  grows 
in  a  multi-rate  architecture. 

The  ratio  between  the  time  frames  i3  two.  Finding  the  denominator 

polynomial  in  the  T/3  time  frame  can  be  considerably  simplified  by  not- 

3  3 

ing  the  denominator  is  a  polynomial  in  z  .  Therefore,  letting  x  =  z 
gives  the  denominator  polynomial,  in  T/6,  as: 

D  »  (x2  -  1. 6x  +  1 ) (x 2  -  1.678134619x  +  .718024891),  T/6 


xl,2  =  .8  +  j. 6  in  T/6  x3  .  +.  839067310  +  j .  1 


18283309 


=*  1 4-  -643501109  x3  4  “  .847363494  4- • 140047156 

2  2 

xj  f  2  =>14-1.28700218  x3>4  -  .  7180248914-  •  280094313 


,28  +  j  •  96 

o’ 


x3  4  “  *690043009  +  j. 198495315 

0 


x2  -  .56x  +  1  in  T/3  x2  -  1.380086018x  +  .515559744  in  T/3 


Thus 


D  -  (z6  -  .56z3  +  l)(z6  -  1.380086018z3  +  .515559744) 

-  z12  -  1.940086018z9  +  2. 288407914z6 

-  1.668799475z3  +  .515559744  (303) 


in  the  T/3  time  frame,  and  also  defines  the  "c"  vector. 


Observe  that 

z6  -  .56z3  +1  -  z3  -  2 (cos  bT)z  +  1  (304) 

represents  the  poles  of  the  sine  wave  input  in  a  T  time  frame  since 
z3  «  (es^/3)3  m  es^. 


To  this  point. 


CT/3 


.014321407 (z9  +  z8  +  z7)  +  .010448513 (z6  +  z5  +  zA) _ 

z12  -  3.278134619z9  +  4.40304021z6  -  2.82697445z3  +  .718024891 


T/3 


(305) 

Dividing  the  numerator  into  the  denominator,  and  saving  every  other 
value  (since  the  ratio  of  time  frames  is  2),  gives  the  "a"  vector  (see 
Table  8) . 

The  answer  will  have  the  form 


cT/3 


12  11 
c^2z  +  cj^z  +•* 


+  Cf 


(z6  -  1.380086018z3  +  .515559744)(z6  -  .56z3  +  l) 


where 


(306) 


c 


12 


c 


11 


o 

CSI 

08 

1 _ 

<S| 

all  3 12 

bn 

all 

•  • 

• 

• 

• 

•  • 

_a0  a0  3 12  _ 

.V 

(307) 
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The  entries  of  the  "a”,  "b"  and  "c"  vectors  are  tabulated  In 

Table  8. 


TABLE  8.  a,  b,  and  c  VECTORS,  PART  III 


a 

b 

C 

a  12 

0.000000000 

b12 

l. 000000000 

c12 

0. 000000000 

all 

0.000000000 

bll 

0.000000000 

C11 

0.000000000 

a10 

0.014321407 

bio 

0.000000000 

c10 

0.014321407 

a9 

0.057396013 

b9 

-1.940086018 

c9 

0.057396013 

a8 

0.057396013 

b8 

0.000000000 

C8 

0.057396013 

a7 

0.125094126 

b7 

0.000000000 

c7 

0.097309364 

a6 

0.197844678 

b6 

2.288407914 

c6 

0.086491475 

a5 

0.197844678 

b5 

0.000000000 

C5 

0.086491475 

a4 

0.24740949 

b4 

0.000000000 

C4 

0.039820806 

a3 

0.259992496 

b3 

-1.668799475 

C3 

0.007502293 

a2 

0.259992496 

b2 

0. 000000000 

c2 

0.007502293 

al 

0.222152094 

bl 

0.000000000 

C1 

0.000000000 

a0 

0.147440916 

b0 

0.515559744 

c0 

0.000000000 

Step  IV; 

T/3 

The  equation  for  C  (Equation  306)  has  now  been  expressed  in  terms 
of  a  sine  wave  Input  in  a  T  time  frame,  since  it  can  be  written  as 
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CT/3 


c12z^  "*■  C1L2^  ***  c0 

sin  bT(z6  -  1.380086018Z3  +  .515559744)z3 


z.esT/3 


_ sin  bTz _ 

z3  -  2  cos  bTz  +1  gT 

z«e81 


(308) 


Following  Section  III  and  Section  VIII,  the  three  components  which  match 
the  output  samples  are  determined  by 


An  +  J®n 


1 


3  sin  bT 


[ci2z^  +  cll2^  +  *' 


+  co] 


[z(z3  -  1.380086018z  +  .515559744)] 


(309) 


z«l  u)nT 


The  precautions  noted  in  Section  VIII  must  be  observed.  The  coeffi¬ 
cients  of  Equation  309  are  tied  to  the  input  frequency  b,  but  the  evalu¬ 
ation  of  z  is  carried  out  in  terms  of  u>n,  where 

-  b  +  ,  n-0,  ±  1  (310) 

To  obtain  a  feel  for  the  effect  of  the  three-rate  sampling  format  on 

the  Bode  plot,  run  a  comparison  with  the  continuous  baseline  case  (Equa- 

T/3  T/3 

tion  286)  by  passing  C  through  a  ZOH.  Thus,  we  may  suppose  the  C 

samples  are  being  reconstructed  into  a  staircase  signal  in  order  to 

drive,  perhaps,  a  motion-base  simulator  actuator.  In  this  event. 


An  +  J®n 


1  -  e-3T/3 
Ts 


s-jwn 


[cl2212  +  ...  +  c0]z_1^t|)nT/3 

sin  bT[z6(z3  -  —)]zml^m 


(311) 


“n 


b 


n 


0,  ±1»  ±2» 


•  •  • 


(312) 
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The  resultant  "primary"  Bode  plot  is  shown  In  Figure  54.  Recall,  from 
Section  V,  that  this  Bode  plot  is  only  one  of  the  two  needed  to  com¬ 
pletely  define  the  coefficients. 

A  most  interesting  phenomenon  is  evident  in  Figure  54,  namely  the 
digitized  system  has  less  phase  lag  than  the  continuous  baseline  system! 
This  is  different  than  the  results  one  would  expect  from  a  single-rate 
simulation  and  therefore  deserves  careful  analysis. 

First  of  all,  the  phase  lead  noted  at  10  rad/sec  cannot  be  attri¬ 
buted  to  a  failure  of  the  computer  program  to  add  in  a  correct  multiple 
of  2ir,  since  we  do  not  expect  to  see  a  phase  difference  of  approximately 
300  deg  at  this  low  frequency  between  the  baseline  and  its  discrete  ver¬ 
sion. 

Next,  a  small  experiment  was  performed  wherein  the  transient 
response  was  run  for  a  sine  wave  input.  The  three  sine  wave  components 
given  by  the  N  =  3  discrete  theory  were  found  to  match  the  steady-state 
output  samples  exactly.  Thus,  the  algorithms  appear  to  be  correctly 
implemented  in  the  computer  code. 

One  might  also  speculate  that  Figure  54  represents  only  one  of  the 
two  Bode  plots  needed  to  completely  specify  the  infinite  set  of  sine 
waves  needed  to  match  the  staircased  signal  that  is  reconstructed  by  the 
ZOH.  Hence,  one  could  argue  that  the  lead  shown  at,  for  example, 
20  rad/sec  might  be  completely  negated  by  the  phase  contributions  of 
other  components. 

As  we  shall  see,  the  same  sort  of  phase  lead  phenomenon  arises  in 
the  larger-scale  case  study,  to  be  discussed  next.  At  this  juncture,  it 
is  felt  that  the  evidence  is  not  strong  enough  to  draw  the  blanket  con¬ 
clusion  that  a  multi-rate  sampling  format  is  an  effective  way  to  intro¬ 
duce  phase  lead  into  a  simulation. 

C.  THE  A-10,  DISPLAYED  PITCH  TO  PILOT  STICK 

FORCE  CASE  STUDY 

This  case  study  represents  an  important  element  extracted  from  an 
overall  simulation  of  the  A-10  aircraft.  It  was  initially  coded  for  use 
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Mognitude 


Figure  54,  A-10  Like  Example 
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Phase  Angle  (deg) 


IP 


in  a  training  (moving-base)  simulator.  Engineers  of  the  Simulator  Sys¬ 
tem  Program  Office  (ASD/ENETS)  proposed  it  as  a  case  study  because  of 
its  interesting  structure.  Specifically,  three  elements  of  the  dis¬ 
played  pitch  angle  to  pilot  stick  force  input  simulation  were  imple¬ 
mented  on  three  different  computers. 

The  initial  rendering  of  the  three-rate  A-10  example  resulted  in  the 
representation  shown  in  Figure  55.  Shown  also  are  the  s-plane  represen¬ 
tations  for  the  eleven  transfer  functions  involved.  The  intent  of  this 
figure  is  to  convey  the  following: 

1.  Each  transfer  function  is  digitized  using  rec¬ 
tangular  integration. 

2.  The  G.  algorithm  feeds  G^,  the  G^  algorithm  feeds 
Gg  and  Gj,  etc. 

The  representation  of  Figure  55  is  a  "shorthand"  which  introduces  no 
difficulties  if  it  is  understood  that  each  transfer  function  is  digi¬ 
tized  individually.  The  Figure  55  representation  would  cause  problems 

if  z-transform  operations  were  to  be  carried  out  since  [g i (s )G2 (s ) ]T  * 
T  T 

Gi(s)  x  G2(s).  To  emphasize  the  individual  digitization  and  the  three 
different  rates,  we  may  redraw  Figure  55  as  shown  in  Figure  56,  where 
data  holds  are  used  to  model  storage  registers.  In  Figure  55,  the  T^ 
data  rate  was  changed  from  156  Hz  to  160  Hz  so  that  a  basic  interval  of 
0.1  second  (as  opposed  to  1.0  second)  could  be  assigned  to  T.  That  is 

1  _L  jJ.  _1  jJ, 

160  “  16  ’  20  “  2  ’  30  3 

so  that,  when  T  »  0.1,  the  rates  can  be  represented  as 

TIT 
11  '  it  •  12  2  '  13  '  3 

At  the  time  this  decision  was  made,  it  was  felt  that  these  rates  were 

sufficiently  separated  to  expose  the  basic  character  of  the  three-rate 

problems  and  to  "simplify"  the  integers  involved. 

Also  indicated  in  Figure  56  is  a  suggested  change  to  the  initial 
drawing.  It  is  probably  realistic  to  model  the  displayed  pitch  angle  as 
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a  staircased  function;  therefore,  a  data  hold,  in  a  T/3  frame  time, 
should  be  inserted  before  the  e  pure  time  delay. 

Since  several  of  the  operations  occur  at  the  same  rate,  we  may  elect 
to  simplify  the  blocks  by  referring  to  Figure  56.  From  the  figure, 
write 


DT/2  -  gi(2[gU2  +  g?/2gV2]ct/2 


(313) 


cl'2 cT/2 


(314) 


T/2  T/2  T/2 
g4  g3  C 


T/2  T/2 
Gb  C 


(315) 


ET/2  -  gI/2gI/2gS/2CT/2  -  GT5/2GTb/2CT/2 


(316) 


AT/2  .  gI/2GT7/2BT/2  -  G5/2BT/2 


(317) 


Using  the  definitions  of  Equations  314  through  317  we  may  draw  the 
more  compact  diagram  of  Figure  57. 

D.  A-10  SWITCH  DECOMPOSITION  HDD EL 


The  switch  decomposition  model  for  Figure  57  is  shown  in  Figure  58. 
We  must  now  inquire  as  to  whether  or  not  the  switch  decomposition  model 
for  Figure  57  would  give  the  same  answer  as  a  switch  decomposition  model 
of  the  more  detailed  Figure  56.  That  is,  the  switch  decomposition  model 
of  Figure  59b  will  give  the  same  answer  as  Figure  59a.  Will  the  switch 
decomposition  model  of  Figure  59c  also  give  the  correct  answer?  The 
answer  to  this  is  in  the  affirmative. 


Figure  57.  4-10  Case  Study-Simplified  Block  Diagram 


Figure  58.  Switch  Decomposition  Block  Diagram 


Figure  59-  Equivalent  Models 


From  Figure  59a 

cT/3  «  [M2G2/2Gi/2] 

and  from  Figure  59b, 

CT/3  -  W3(w3*M2W2)T(w2*G2W2)T 

x  (w2*g1w2)t(w2*r)t 

From  Figure  59c,  one  finds 

CT/3  -  W3(w3*M2W2)T(W2*GW2)T(W2*R)T 

Now,  rearrange  Equation  320  (nest  the  T  operators) 

CT/3  -  {[w3(w3AM2W2)(w2(tGW2)(w2^R)T]T} 

Notice  W2(W2AR)T  -  RT/2,  therefore 

cT/3  -  {[w3(w3*m2w2)(w2agrt/2]t} 


Next,  note  W2(W2*GR^/2 is  by  definition  GT^RT^  which  is,  in  turn 
G2^2G^2R^2.  The  last  T  operator, 


T /3 

cT/3  „  W3(w3j|[M2G2/2Gi/2Rt/2)T  -  [m2G2/2g]/2RT/2  ]  (322) 

Thus  the  switch  decomposition  model  of  Figure  59c  is  exactly  equiva¬ 
lent  to  the  model  of  Figure  59b. 

When  this  study  was  initiated,  switch  decomposition  appeared  to  be 
the  logical  solution  to  the  problem.  However,  the  developments  of  Sec¬ 
tion  III  pushed  the  decomposition  approach  into  the  background.  Fig¬ 
ures  58  and  59  are  therefore  given  only  to  illustrate  an  alternate 
approach. 

E.  ALGEBRAIC  MANIPULATIONS 

Referring  to  Figure  57,  write 


AT/2  =  -  [M!Ga/16RT/16]T/2  +  {HiGa/16[m2Gc/2At/2] 


T/16 


T/2 


♦  Gl/2GpAT/2  (323) 


Simplify  Equation  323: 

T/2  T/2 

aT/2  .  _  [MiGA^I6RT^16]  +  iH1CTA/'V/16  ]  g£/2At« 

♦  g£'2g5/2aT/z  (324) 

Therefore, 

A1"  -  -  (I  -  [KlGi[/16Hl/16lT/2Gj/2  -  g5/2gE/2!_1 

T/2 

*  [MlGA/16RT/i6]  (325) 
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The  final  result  is 


T/3 

0  =  G6|M2G5/2Gb/2At/2]  (326) 

or 

0  =  -G6{m2G5/2Gb/2[I  ‘  (MiGa/16M2/16)  '  -  g5/2G^/2] 

T/2 

x  [MiGj^ 16rT/16 ]  }  (327) 

Comparing  this  with  Equation  277,  it  can  be  appreciated  that  the 
A-10  case  study  reduces  to  the  same  type  of  computations  required  for 
the  simpler  example.  The  dimensionality  will  be,  of  course,  consider¬ 
ably  higher.  We  will,  therefore,  not  repeat  the  analysis  of  the  pre¬ 
vious  section,  but  merely  note  that  a  computer  program  was  written  for 
Equation  327  and  the  results  computed  on  a  Tymshare  terminal.  The  s- 
and  z-plane  transfer  functions  are  listed  in  Table  14  in  Appendix  D.  We 
proceed  immediately  to  a  discussion  of  the  results.  To  compare  the 
difference,  the  continuous  baseline  Bode  plot  is  shown  with  the  Bode 
plot  obtained  by  assuming  that  the  displayed  pitch  angle  sample  sequence 
is  reconstructed  with  a  ZOH  (in  order  to  drive  the  display  generator). 

F.  COMPARISON  WITH  BASELINE 

The  Bode  plot  associated  with  the  three-rate  simulation,  together 
with  that  of  the  continuous  baseline  system,  is  shown  in  Figure  60. 
Again,  one  must  keep  in  mind  that  this  is  only  one  of  the  two  plots 
needed  to  completely  define  the  spectral  components  (although  it  does 
appear  to  be  the  primary  one). 

As  was  the  case  with  the  lower-order  examples  used  previously,  we 
note  the  phase  lead,  at  low  frequencies,  is  less  for  the  digitized  ver¬ 
sion  than  it  is  for  the  continuous  model.  Clearly,  this  phenomenon 
requires  further  study  and  we  caution  against  adopting  the  viewpoint 
that  a  multi-rate  sampling  format  is  a  technique  which  will  always 
introduce  phase  lead  into  a  simulation. 
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Frequency  trad/sec) 


Figure  60.  Display  Pitch  Angle  to  Stick  Force  Input 


C.  SUMMARY  AND  CONCLUSIONS 


Two  case  studies,  which  utilized  three  different  sampling  rates, 
were  analyzed.  Each  introduced  phase  lead  into  the  first  fold  of  the 
frequency  response,  a  result  which  contradicts  our  intuitive  expecta¬ 
tions,  It  was  shown  that  the  multi-rate  analysis  tool,  developed  in 
Section  IV,  was  an  effective  tool  for  treating  the  dimensional  problems 
Introduced  by  the  variations  in  the  three  sampling  frequencies. 


i 

! 

I 


SECTION  X 


AN  ORDER-OF -CALL  CASE  STUDY 


A.  INTRODUCTION 

This  section  focuses  on  a  single-rate  case  study  which  is  concerned 
with  order-of-call  effects*  These  anomalies  are  introduced  when  one 
element  of  code  calls  another  element  containing  results  predicated  on 
"fresh"  input  data  but  which,  in  realitv.  used  "stale"  data  (because  the 
new  data  were  not  yet  available).  Specifically,  we  study  a  roll/sway 
washout  network  simulation  wherein  a  pure  delay  is  introduced  in  an 
integration  algorithm.  This  case  study  was  suggested  by  AFWAL/FIGD  per¬ 
sonnel  who  are  associated  with  LAMARS.  The  problem  was  first  described 
in  Reference  5.  The  following  description  of  the  problem  is  excerpted 
from  that  reference. 

B.  DESCRIPTION  OF  WASHOUT  COMPUTATION 

In  the  LAMARS  motion  washout  routine,  transfer  functions  are  imple¬ 
mented  by  representing  the  output  variable  as  a  function  of  the  input 
variable  and  also  representing  it  implicitly  as  integrals  of  the  output 
variable  itself.  The  output  equation  is  evaluated  first,  then  the 
derivatives  of  the  states  are  computed  based  upon  the  output  and  lower 
states,  then  finally  all  integrations  are  performed  using  the  Adaros- 
Bashforth  second-order  predictor  algorithm. 

For  a  third-order  washout  filter  of  the  form  used  in  the  w  path,  the 
implementation  is  developed  as  follows: 


2 

Output  _ fa  3s _ 

Input  r  s^  +  a2S^  +  ajs^  +  a^ 


(328) 
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3 

[Divide  through  by  s  ] 


1 

r 

[Cross  multiply] 

a  a  a 

(l  +  -2.  +  -i  +  -f)  y  -  b3r  (330) 

S  gZ  gj 

[Move  integrations  to  input  side] 

y  -  b3r  -  — •  [a2  +  ~  (ai  +  ~  a^J  ]y  (331) 

[Define  output  in  terms  of  input  and  highest-order  state] 

y  *»  b3r  -  w3 

[Define  states  from  deepest  nesting  outward] 

wi  -  a^y  wi  -  -  Wi 

w3  “  aiy  +  ~  a<j>y  “  aiy  +  wj  — ^  w3  «  ~  W2  (332) 

s  T  s 

w3  -  a2y  +  —  aiy  +  —  2a*y  -  a2y  +  W2  w3  -  —  w3 

Thus,  the  input  is  seen  immediately  in  the  output  and  thence  in  all 
derivatives  of  the  states,  but  the  states  themselves  are  not  updated. 
The  expansion  outward  from  the  nesting  results  in  state  equations  which 

assume  in  their  form  that  updated  states  are  continuously  available. 

However,  since  all  loops  are  then  implicitly  expanded  in  full  prior  to 
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any  integrations,  a  lag  with  respect  to  the  system  state  response  is 
introduced,  while  immediacy  of  response  with  respect  to  the  input  is 
preserved. 

C.  THIRD-ORDER  WASHOUT  CASE  STUDY 

The  analog  block  diagram  of  the  previous  equations  is  given  in  Fig¬ 
ure  61b  and  then  translated  into  the  digital  version  of  Figure  61a.  In 

Figure  61a,  there  is  an  explicit  model  of  the  order  of  call.  That  is, 

the  integrations  are  performed  last  in  the  cycle.  This  has  been  inter¬ 
preted  to  mean,  on  any  given  machine  cycle,  that  wj>  w2>  ate  the  old 

values  (previous  frame  time  values).  Thus,  the  integration  of  w^  that 

contributes  to  w2  is  one  frame  time  old,  the  integration  of  w2  which 

contributes  to  W3  is  a  frame  time  old,  as  is  the  integration  of  W3  which 

produces  w^> 

In  evaluating  (l/s)T,  Reference  5  calls  for  the  use  of  an  Adams- 

Bashforth  second-order  predictor  algorithm.  It  was  later  agreed  that 

nothing  essential  would  be  lost  by  using  "Implicit"  Adams-Bashforth 
second-order,  which  is  described  by  Figure  62  and  the  equations 


Xn_!  +  (T/2)[3Xn  -  Xn_i] 


(333) 


which  gives. 


z  -  I 


(334) 


If  the  integrations  could  be  accomplished  with  very  little  throughput 
delay,  T0  -  0,  one  could  write,  after  setting 


z  -  1 


-  b3RT  -  [a2G  +  aiG2  +  aQG3]yT 


(335) 


(336) 


i  _ .. 

,  .  .. , 


•*  .#■  .** 


Figure  62.  Discrete  Integrator 


or 


[l  +  a2G  +  a^G2  +  a0G3]-1b3RT 


(337) 


With  time  delay,  the  equation  for  W2  in  terms  of  W2  becomes 

-sT  o 


W2 


(338) 


If  Tq  -  T,  this  reduces  to 


W2 


(T/2)  (3z  -  1)  ‘ 


z(z  -  I) 


W2 


(339) 


G  is  modified  by  the  factor  z-*.  This  procedure  can  be  repeated  on  the 
other  two  "nestings"  giving  the  same  result.  Each  G  of  the  no-delay 
case  is  only  modified  by  z-1  and  Equation  347  can  be  used  for  both 
cases. 

The  result  obtained  requires  careful  validation  since  it  forces  the 
characteristic  equation  of  Equation  337  to  increase  by  an  order  of  two. 
To  see  this,  insert  the  appropriate  G(z)'s  and  clear  through: 


{[i  +  (3/2)Ta2  +  (9/4)T28i  +  (27/8)T3a0]z3 

+  [-3  -  (7/2)Ta2  -  (I5/4)T2ai  -  (27/8)T3a0]z2 

+  [3  +  (5/2)Ta2  +  (7/4)123!  +  (9/8)T3aQ]z 

+  [-1  -  (T/2)a2  -  (T/4)2ai  -  (T/8)3a0]}  yT  -  b3(z  -  1)3rT 

No  Delay 


(340) 
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{z6  +  [ -3  +  (3/2)Ta2]z5 


+  [3  -  (7/2)Ta2  +  (9/4)t231]z4 
+  [-1  +  (5/2)Ta2  -  (15/4)T2ai  +  (27/8)T3a0]z3 
+  [-(7/2)a2  +  (7/4)t23i  -  (27/8)T3a0]z2 

+  [-(7/4)2ai  +  (9/8)T3ao]z  -  (7/8)3a()}  “  b3z3(z-l)3RT 

(341) 

Delay 

For  G  =  (6s^)/(s^  +  6s^  +  11s  +  6)  and  T  “  0*04,  Equation  340  yields 
three  real  roots,  all  interior  to  the  unit  circle  (we  are  using  conven¬ 
ient  numbers;  more  realistic  parameter  values  will  be  used  for  a  sixth- 
order  case,  which  follows). 


A(z)  *  (z  -  .898304953) (z  -  . 962264027) (z  -  .928571684)  (342) 

With  the  frame  time  delay  in  the  integration.  Equation  351  gives 
A(Z)  -  (z  -  . 887598) (z  -  . 960815504) (z  -  .92332214) 

(343) 

x  (z  +  .067598089) (z  +  .043321992) (z  +  .020815562) 

The  results  for  the  sixth-  and  third-order  models  are  summarized  in 
Equations  344a  and  344b,  the  data  for  which  are  given  in  Table  9. 


Gl 


ap (z3) (z  -  l)3 


+  b5a3  +  b^z4  +  b3Z3  +  b2z2  +  bjz  +  b( 


(344) 


or 


Cn(z  -  1) 

Q «  ■  ——■■■■  '  ■■  ■■■ 

z^  +  d2z2  +•  d^z  +  dg 


(345) 
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TABLE  9.  THIRD,  SIXTH  MODEL  DATA 


bo 

-0.000048000 

d0 

-0.802663438 

bl 

-0.003968000 

dl 

2.592078213 

b2 

-0.090496000 

d2 

-2. 789140664 

b3 

-0.464704000 

c0 

4.282973183 

b4 

2. 199600000 

*>5 

-2.640000000 

ao 

6.000000000 

The  results  of  this  analysis  are  most  interesting.  They  indicate 
the  order  of  call  has  introduced  an  effective  filter  which  essentially 
doubles  the  order  of  the  analytical  model  describing  the  software.  In 
the  illustrative  example  we  see  that  the  "delay"  case  tracks  very  well 
with  the  continuous  baseline  and  no  delay  case  (See  Table  10).  Indeed, 
one  could  conclude  that  the  delay  case  is  superior  to  the  no  delay  case. 

The  correctness  of  the  result  can  also  be  checked  by  running  the 
impulse  response  using  both  the  analytical  model  of  the  simulation  and 
the  actual  computer  code.  The  results  shown  in  Table  11  indicate  an 
extremely  close  agreement. 

Thus  we  see  that  the  manner  in  which  an  algorithm  is  implemented  can 
introduce  filtering  over  and  above  that  intended  by  the  designer.  We 
shall  explore  this  in  more  depth  by  next  considering  a  sixth-order  wash¬ 
out  using  parameter  values  supplied  by  FIGD  personnel. 

D.  SIXTH-ORDER  WASHOUT  EXAMPLE 

Here  we  treat  a  sixth-order  model  where  the  transfer  function  repre¬ 
sents  lateral  acceleration  (y^)  to  roll  rate  (PA): 
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TABLE  10.  SOME  REPRESENTATIVE  BODE  DATA 


TABLE  11.  IMPULSE  RESPONSES 


k 

6th  ORDER 
MODEL 

COMPUTER 

CODE 

6.000000000 

6.000000000 

1 

3.840000000 

3.840000000 

2 

2.940000000 

2.940000000 

3 

2. 103360000 

2. 103360000 

4 

1.413485760 

1.413485760 

5 

0.842594150 

0.842594150 

6 

0.374368440 

0.374368444 

7 

-0.005989020 

-0.005989006 

8 

-0. 311312900 

-0.311312875 

9 

-0.552761020 

-0. 552760985 

10 

-0. 740018250 

-0. 740018208 

11 

-0.881479430 

-0.881479389 

12 

-0.984410188 

-0.984410140 

13 

-1.055088400 

-1.055088363 

14 

-1.098928750 

-1.098928718 

15 

-1.  120592090 

-1. 120592066 

16 

-1.  124081650 

-1.124081640 

17 

-1.112827510 

-1. 112827520 

18 

-1.089760760 

-1.089760787 

19 

-1.057378550 

-1.057378593 

20 

-1.017801150 

-1.017801208 

21 

-0.972821940 

-0.972822016 

22 

-0.923951200 

-0.923951295 

23 

-0.872454410 

-0.872454524 

24 

-0.819385750 

-0.819385882 

25 

-0.765617360 

-0. 765617512 

26 

-0. 711864890 

-0. 711865060 

27 

-0.658709770 

-0.658709954 

28 

-0.606618610 

-0.606618808 

29 

-0.555960090 

-0.555968303 

t  -  k  At,  At  «=  0.04  sec 


■> 
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1.61s’  8  +  1.56923 


[s2  +  .07s  +  .0025  ] (s  +  1.57)(s3  +  1.56923s2  +  s  +  .266) 


1.61s  +  2.5264603s* 


r30s  +  R298  +  R288  +  R278  +  R268  +  R258  +  R24 

where  the  R^j's  are  listed  In  Figure  63.  We  may  nest  Equation  347  as 
follows: 


The  analog  block  diagram  is  shown  in  Figure  63. 

We  may  dispense  with  the  digital  block  diagram  and  write  the  appro 
priate  equations  directly  from  Figure  63. 


1)  Y  -  Wo  -  W 


Wi 

rx 

r24Y 

3)  Wi 

-  T/2  — 

z  — 

w2 

“ 

*25* 

+ 

Wi 

i  «  1,  •  • 

w3 

m 

R26* 

+ 

W2 

• 

W4 

m 

R27Y 

+ 

W3 

or 

w5 

m 

r28y 

+ 

W4 

w, 

ik 

-  W. 
ik-1 

W5 

m 

r28y 

+ 

W4 

+  3T  /  2Wij 

«6 

m 

r29Y 

+ 

W5 

-  T/2Wik 

«7 

m 

R£2r 

i-1,  2, 

W8 

m 

r23r 

+ 

W7 

'•  «r* 


Figure  63.  A  More  Complex  "Sway/Roll"  Example 


G(z),  for  the  no-delay  case  Is  readily  computed  and  has  the  form 


x 

R 


R11926  +  Rll8z5  +  Rj 17Z4  +  r116z3  +  Rll5z2  +  R114z  +  Rl 1 3 
r12626  +  r125z5  +  r12424  +  R123z3  +  R222z2  +  R221z  +  R220 


and  for  the  delay  case 

Y  _  rI12212  +  RlllzU  +  r11Qz10  +  +  Rj00  (352 

R  r22z^2  +  r21z^  +  r20z^  +  ***  +  R10 

The  data  are  tabulated  in  Table  12  (again,  T  «  0.04  seconds). 

Again,  we  found  that  the  impulse  responses  of  the  twelfth-order 
analytical  model  agree  exactly  with  the  impulse  response  obtained  dir¬ 
ectly  from  the  computer  code  (not  shown). 

Some  representative  Bode  plot  data  are  given  in  Table  13. 

Using  the  parameter  values  given  us,  we  see  that  the  filtering 
Introduced  by  the  order  of  call  has  actually  improved  the  fidelity  of 
the  frequency  response. 

B.  SUMMARY  AND  CONCLUSIONS 

The  encoding  of  an  algorithm,  where  a  definite  sequence  of  call  is 
given,  can  introduce  filtering  of  which  the  designer  of  the  algorithm  is 
unaware.  For  the  washout  case  study,  we  see  that  the  filtering  action 
was  "benevolent"  in  that  it  actually  improved  the  fidelity  of  the  fre¬ 
quency  response.  The  effect  of  updating  the  integrator  states  as  the 
last  event  in  a  given  frame  time  leads  to  an  analytical  model  of  the 
computer  code  which  is  twice  the  order  of  algorithm  being  encoded. 
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TABLE  12.  COEFFICIENTS  FOR  6th  ORDER  WASHOUT 
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TABLE  13.  FREQUENCY  RESPONSE  DATA,  6th  ORDER  WASHOUT 


MAG 

i- 

freq. 


CONT 

DISCRETE, 

NO  DELAY 

DISCRETE 

DELAY 

1.833472489 

-2.614690348 

1.948711343 

4.026912923 

36.32463432 

14.36929036 

1.000000000 

1.000000000 

1.000000000 

-2.405029014 

-3.026181390 

-2.412522944 

-42.92847068 

-38.25338109 

-42.44121200 

2.000000000 

2.000000000 

2.000000000 

-5.633216300 

-6. 144060048 

-5.576132340 

-58.65023633 

-52. 10560145 

301.5580652 

3.000000000 

3.000000000 

3.000000000 

-8.031772688 

-8.504798812 

-7.938527142 

-66.49678734 

-57.69456014 

293. 7055865 

4.000000000 

4. 000000000 

4.000000000 

-9.924000906 

-10.34754451 

-9.778310534 

-71.20086176 

-60. 12882875 

289.0036661 

5.000000000 

5.000000000 

5.000000000 

-11.48275868 

-11.84540095 

-11.27354613 

-74.33561069 

298.9727785 

-74. 15427992 

6.000000000 

6. 000000000 

6.000000000 

-12.80673528 

-13.09802451 

-12.52335630 

-76.57418687 

298.9268998 

-76.44184967 

7.00000000 

7.000000000 

7.000000000 

-13.95687896 

-14.  16724547 

-13.58926581 

-78.25287151 

299.4021514 

-78. 20508987 

8.000000000 

8.000000000 

8.000000000 

-14.97328861 

-15.09375872 

-14.51203844 

-79.55838679 

300.2095026 

-79.63326613 

9. 000000000 

9.000000000 

9.000000000 

15.88369136 

-80.60272779 

10.00000000 


-15. 90583247 
301. 2347408 
10.00000000 


-15.31989859 

-80.84148715 

10.00000000 


APPENDIX  A 


EXAMPLES  DEMONSTRATING  THE  SETUP  OF  THE  GENERALIZED 

SKIP-SAMPLING  THEOREM 


This  appendix  contains  several  exercises  on  the  setup  procedure  for 
the  generalized  skip-sampling  theorem.  The  situation  is  reviewed  in 
Fig.  A-l  and  Eqs.  A-l  and  A-2. 

CT/M  „  [GrT/N]T/M  =  [GT/NKRT/N]T/M  (A-i) 

In  Eq.  A-l,  M  and  N  are  rational  numbers  and  K  is  an  integer.  If 
z  =*  esT/NK^  then  NK/M  must  be  an  integer.  Equation  A-l  can  be  evaluated 
using  Eq.  A-2: 

cT/M  -  7FJ  jf  <=(P)«PK)  z  rjrnn  %  ■  » - 

(A-2) 

V"*  . ,  ,  G(p)R(pK) 

=  3’  residues  of  - - 


Figure  A-l,  A  Phantom  Sample  Formulation  of  a 
T/N,  T/M  Sampling  Format 
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APPENDIX  B 


SKIP-SAMPLING  THEOREM 


When  dealing  with  functions  such  as 


[grt/n]T/M  =  [gt/mmrt/N]T/m 


(B-l ) 


there  is  the  option  of  which  poles  to  circle.  That  is,  in 


„  T/M  ,  r 

[GRT/N]  =  jij  Jr  G(P>R(Pm> 


z  dp 
z  -  pn  P 


(B-2) 


N  N 

the  poles  of  G(p)R(p  )/p  or  the  poles  of  z/(z  -  p  )  may  be  enclosed. 
For  the  special  case  of  [rt/n]^,  the  result  of  closing  the  contour 
around  z/(z  -  pN)  is  known  as  Sklansky's  identity  (Ref.  6).  The  result, 
for  arbitrary  N,  is 


[*T/fT  -  i  E  *(%«■) 

k-0 


where 


zN  -  eaT/N 


ak  -  eJ2nk/N 


(B-3) 


k  -  0,  1,  ...,  N  -  1 


Of  particular  interest  is  the  case  of  N  ■  2,  since  only  real  numbers 
are  involved: 


afc  - 


1  k-0 


k-1 


1  k-0 


-1  k-1 


[RT/2J  -  ■}■  [R(2>  + 


z  —  esT^ 


;  V  V  '  L&Jit 


(B-4) 


The  usefulness  of  Sklansky's  identity  for  frame  time  ratios  which  are 
powers  of  2  is  clear.  For  example,  if 


CT  .  [rT/16]T  . 

f[RT/16]T/S 

L 

T/4 

T/2 

- 

This  problem  can  be  solved,  using  only  real  numbers,  by  applying  Sklan¬ 
sky's  Identity  four  times  (2^  ■  16). 
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APPENDIX  C 


A  DIMENSIONALITY  PROBLEM 


The  use  of  the  high-to-low  rate  transfer  often  introduces  a  dimen¬ 
sionality  problem.  For  example,  if 


cT/3 


M^6  GT/2 


T/3 


T/3 


z2  ±  z 


±1 


G  (z  3) 


(C-l) 


zme 


sT/6 


3 

and  G(z)  is,  for  example,  a  seventh-order  transfer  function,  G(z  )  is 
21st-order  and  therefore 


[•] 


z2  +  z  +  i 
Z2 


G  (z  3) 


is  23rd*  This  dimensionality  is  of  little  concern  if  one  chooses  to 
implement  the  algorithm  of  Section  IV;  however,  the  occurrence  of  shift¬ 
ing  between  time  frames  ratioed  by  powers  of  2  is  common  enough  to  make 
a  technique  for  using  Sklansky's  identity,  in  a  more  efficient  manner, 
attractive. 

Suppose  we  replace  z^  by  x  in  G(z^)  and  find  the  high-to-low  rate 
transform,  with  a  ratio  of  2.  That  is, 

[G(x)]t/2  -  W(x) 


i 

t 


\ 


Using  Sklansky  s  identity,  we  may  write: 


1  I  z2  +  z  +  1  _  .  -k  .  z2  -  z  +  1 
_  - - -  G(zJ)  + - r - 


1  G(z 


2-i-iW(z6)  +  7[ 


G(z  )  -  G(-z ' 
2 


(C-4) 


The  reader  may  find  it  puzzling  that  a  z  argument  is  placed  on  W.  To 
see  this,  recognize  that  the  use  of  Sklansky's  identity, 


W<22>|z=esT/2  -  w(2>|z=esT 


(C-5 ) 


generates  a  function  of  z  .  That  is,  when  one  evaluates  equations  like 

2  sT  /  2 

Eq.  C-5,  we  are  dealing  with  a  function  of  z  (z  =  e  )  and  it  is  con- 

sT 

venient  to  switch  directly  over  to  a  function  of  z  wherein  z  =  e 

However,  in  Eq.  C-5,  recognize  that  the  algorithm  automatically  makes 

2 

the  substitution  for  us.  That  is,  it  does  not  yield  W(x  )  in  the  T/2 

2  2 

time  frame;  rather  it  gives  W(x)  in  a  T  time  frame.  Since  (z  +  1  )/z 

is  still  in  a  T/6  time  frame,  one  must  insure  that  the  arguments  are 

2  3 

compatible.  Clearly,  W(x  )  is  the  desired  format.  Since  x  =  z  ,  we 

have  W(z6). 

Next,  suppose 


H(z3) 


(C-6) 


C-2 


and  apply  the  definitions  to  find  out  what  H  Is,  In  terras  of  G  and  W. 
Since 


G  -  G 
2 

G  +  G 


the  sum  gives 


H  +  W 


(C-7 ) 


Using  Eq.  C-7,  rewrite  Eq.  C-4  as 


z“_+z  +  l 


G  (z3 ) 


T/3 


(C-8 ) 


where 


«  W 


The  use  of  Eq.  C-8  is  clear  —  we  may  now  use  any  convenient  compu¬ 
ter  multiply  option  and  add  option  to  generate  the  23rd-order  CTy/3  with 
a  minimum  of  dimensionality. 

Alternatively,  for  frequency  response  purposes,  we  can  evaluate 
CT/3,  for  z  *  l4bT/6,  rather  easily: 

c^U-ubT/e  -  ~  ■ 


z-4bT/6 


w(z> |z«l*bT 


+ 


I 

z 


z-UbT/6 


G(Z)  lz-UbT/2 


(C-9) 


APPENDIX  D 


s-  AND  z -DOMAIN  A-10  TRANSFER  FUNCTIONS 


The  s-  and  z-plane  A-IO  transfer  functions  used  in  the  A-10  case 
study  are  tabulated  in  Table  14.  The  appropriate  time  frame  used  in  the 
discretization  of  each  transfer  function  is  listed  in  the  left-most 
column  of  the  table. 
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